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An Investigation of the Detached Shock in 
Front of a Body of Revolution’ 


JOHN DUGUNDJIt 


Grumman Aircraft Engineering Corporation 


SUMMARY 


Some theoretical investigations are made into the nature of the 
flow behind the detached shock formed in front of a body of revo- 
lution. The predicted values of the product of the curvature of 
the shock and distance of the shock in front of the body agree 
reasonably well with the experimental values of Charters and 
those of Ladenburg, Van Voorhis, and Winkler. 


(1) INTRODUCTION 


5 ee INVESTIGATION carried out in this article is a 
continuation of that by Lin and Rubinow! and 
deals specifically with the problem of the detached 
shock formed in front of a body of revolution placed 
ina uniform supersonic stream. As has been indicated 
by Lin and Rubinow, there may, in general, be ex- 
pected to be some approximate interrelation between 
the following three parameters: (a) curvature of the 
shock wave at the nose of the shock, 0a/0e; (b) dis- 
tance from the shock to the body along the axis of sym- 
metry, 6; (c) curvature of the nose of the body, ko. 
These relations are now investigated in some detail 
and compared with some experimental results. 

The general plan of such an investigation has been 
discussed in reference 1 and is briefly sketched here. 
One calculates the spatial derivatives of pressure or ve- 
locity along the axis of symmetry in terms of the curva- 
ture of the nose of the shock (§ 2). A power series ex- 
pansion can then be made. If the distance between the 
shock and the body is small, one may use this series to 
calculate the flow conditions at the nose of the body. 
By using known relations for the conditions there, the 
desired interrelation can be found in the form, 


Received May 26, 1948. 

* The author is indebted to Prof. C. C. Lin for suggesting the 
problem and for his guidance and supervision of the work. He is 
also indebted to Prof. H. S. Tsien for his suggestion of the ap- 
proach used in (§2). 

Aeronautical Engineer. 


5(Qa/de) = kod) 


In the above, MM, represents the free stream Mach 
Number. The numerical results thus obtained agree 
satisfactorily with experimental observations (§ 3). 
These and some related quantities are presented in Figs. 
1-5 later in this paper. 

The separate determination of 6 and 0a/0e can be 
made by the additional use of the power series expan- 
sion starting from the nose of the body, as was also out- 
lined in reference 1. This is discussed in (§ 4). How- 
ever, the results are not completely satisfactory, since 
the accuracy is limited by the number of terms in the 
series available at the present. 


(2) DERIVATIVES OF PRESSURE BEHIND THE SHOCK 


The purpose of this section is to obtain some relation 
between the spatial derivatives of pressure immediately 
behind the shock and the curvature of the nose of the 
shock. In the development that follows, subscript 2 
will refer to conditions immediately behind the shock, 
subscript 1 will refer to conditions immediately in 
front of the shock, and subscript 0 will refer to condi- 
tions at the nose of the body. 


(a) First Derivative 


The basic equations for our discussion are the two 
dynamic equations, the continuity equation, and the 
equation of state. 


ou Ou 10p 
1 
pox (1) 
ov ov 10p 
(2) 
re) 
5 pu) (pv) € (3) 


CAS 
CoR- : 
| 
\ 
ON 
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Op/Op = = yp/p (4) 


In the above, u and v represent the velocities in the 
x and y directions, respectively; p, the density; ), 
the pressure; c, the acoustic velocity; and y, the ratio 
of specific heats (y = 1.4). The ¢ of Eq. (3) is a con- 
stant used to make the above equations applicable to 
both the two-dimensional case (« = 0) and the axially 
symmetric case (€ = 1). 

Along the axis of symmetry, 


v= oF 5) 

dv/dx = 0 

By applying the basic equations to the region im- 

mediately behind the shock and substituting condi- 

tions (5) into Eqs. (1), (2), and (3), Eq. (3) can be de- 
veloped into an expression of the form, 


Ox 2 oy 2 


This equation expresses the pressure gradient in 
terms of (Ov/Oy)2. To find the value of (Ov/Oy).s, the 
shock relations will be used. 

The vertical component of velocity behind the shock 
can be expressed as 


= sin (7) 


In the above, g represents the velocity along a stream- 
line, and @ represents the turning angle through the 
shock. 

If Eq. (7) is differentiated along the shock (¢) and 
the condition of normal shock (@ = 0°) is applied, one 
obtains 


(Ov/Oc)2 = 00/00 (8) 


To obtain 00/0c, two shock relationships are intro- 
duced—-namely, 


gi COS @ = Gq. cos (a — A) (9) 


{2 + (y 1) M,? sin? al/(y 1) 
sin? a (10) 


In the above equations, a represents the angle be- 
tween the approaching streamline and the shock. 

Eq. (9) is differentiated with respect to o, the normal 
shock conditions 6 = 0°, a = 90° are introduced, and 
the resulting equation is solved for 00/Oe. This 00/00 
is then substituted into Eq. (8). With the aid of Eq. 
(10) and the fact that 0v/Oo — dv/dy at the axis of 
symmetry, the value of (0v/Ov). becomes 


2 


1) Me de 
By merely substituting this into Eq. (6) and then 
making use of the continuity relation, 


Pilly = pole (12) 


one arrives at the final result, 


(13) 


(M,? — 1) (=) _ 2(1 + (Mi? — 1)0a 

pill,” Ox/ +1 Oa 

Eq. (13) with e = 0 is seen to be in agreement with 

that derived by Lin and Rubinow bv using intrinsic 
coordinates. * 


(b) Second Derivative 

To arrive at the relations for the second derivative 
of pressure, one must differentiate the starting rela- 
tions, Eqs. (1), (2), and (3), with respect to x, y, and 
x, respectively, and then introduce conditions (5), 
The differentiated continuity relation will then be 
changed to 


(M,? — 4 + 


polls? Ox 


(14) 


If one works on the last term of the above equation by 
using Eqs. (11), (13), and (10), the last term will be- 
come 


Op 
(1+ = (=) (22) = 


(y + 1)M.? (1 — 1 (22) (15) 
(2+ 1)M,?] pitty” \Ox/» 


The only remaining difficulty is the (0°v/Oxdy)2 term. 
This can be related to the ordinary variables through 
the use of the differentiated Eq. (2). In this equation, 
(0°v/Oxdy)2 can be solved for in terms of (Ov/Oy)s and 
(0°p/Oy?)2, both of which can be determined from the 
shock relations. The value of (Ov/Oy)2 has already been 
presented in Eq. (11), while the value of (0°p/Oy*). can 
be obtained by differentiating the pressure gradient 
along the shock, then setting a = 90°, thus, 


dp 0a Op p Op 
= cos @ ax + sin a ay 
da? (22) dade? do? 
oO op OP 
(cos a + sin )(cos a + sina 
Oy*/2 Ox/ 2 Oe 


In the above, the 0°a/0c* term drops out, since, by 
symmetry considerations, 


* See Eq. (3.12) of their paper, with 6; = 0°, a = 90°. Use 
Eq. (10) of the present paper to carry out of reduction of g2/q- 
The method used here was suggested by Prof. H. S. Tsien. It 
is more convenient for handling the conditions of the nose of the 
shock. 


Also, 
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= = = 0 


From the shock relations it is known that, 
1 
M,? sin? a + (17) 


(= 


Eq. (17) is differentiated twice, and the result is sub- 
stituted into Eq. (16). The value of (0°v/Oxdy), is now 
completely determined in terms of useful variables. 
Substituting these values into Eq. (14) gives the follow- 


ing basic equation: 
+ {y — 


(M2? — 1) - 
pill,” 
(1 
(22) +(1+6)x 


1 | 
(2 2 ox 
4 ee)’ 4 


+ (y — 1) \o 


(y 1)M,° *(y M2! + 1) 
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Eq. (18) with e = 0 is seen to be identical in the 2- 
dimensional case with Rubinow’s results.* 


The higher derivatives can be determined in much the 
same manner as that described above. One merely has 
to differentiate the basic starting relations [Eqs. (1), 
(2), and (3)] two or more times and then proceed as 
before. 


It is more convenient to introduce dimensionless 


quantities. Thus, Eqs. (13) and (18) can be rewritten 
in the form, 
(1/po)(Op/Ox)2 = u(Oa/do) (19) 
(1/po)(0*p/Ox")2 = v(Oa/da)* (20) 
where, for the axially symmetric case, 
* Eq. (38) of Rubinow’s Thesis. 
(21) 


Sy (1 + [(y — — + 1) + (1 — 


— 1)? [2yM? — (y — 1)] 


Also, the quantity 7 is introduced. 


(pe Do) / Po 


(3) SERIES EXPANSION AROUND THE NOSE OF THE 


SHOCK 


Either the velocity or the pressure can be expanded 
as a power series around the nose of the shock. It was 
found, however, to be more convenient to approximate 
the velocity rather than the pressure distribution prin- 
cipally because of its nonvanishing slope at both the 
shock and the nose of the body.t 

The velocity at any point along the axis of symmetry 
behind the shock can be expressed approximately by 
the power series 


Ou 1 
mt + + 


1/o*u\ 
ax? (24) 


In the above, u is the velocity a distance x from the 
Shock, and (0u/Ox)s2, (O07 /Ox*)o, and (0%u/Ox*)s are 
the velocity and velocity gradients immediately behind 
the shock. At the nose of the body, located a dis- 


k= 


+ A more detailed discussion may be found in reference 2. 


{1 + [(y — 1)/2] 1 


(M2 — 1)[2 + — 
(M2 — 1) + M2 (1 (22) 


(23) 


tance x = 6 behind the shock, the velocity is equal to 0. 
Therefore, Eq. (24) applied to the nose of the body, 
after having been multiplied through by pet%2/Po, and, 
after having had its velocity gradients related to the 
pressure gradients through the Bernoulli equation, 


Op/Ox = —pu(Ou/Ox) (25) 
will give the following expression, 
Oo Oa 
where 
g = yM2*(r + 1) (27) 
= — + — (28) 


Figs. 1 and 2 give values of g, h, and wu for a Mach 
Number range of 1.25 to 4.50. 

If one uses two terms of the series, Eq. (26), one will 
obtain the following expression for 6(0a/0«): 
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—) = 
Po dx 0 Oa 
52 5 1 
(2) b+... (31) 
44 4 
It only remains to determine (du/dx)o, substitute it 
a Ls into Eq. (31), and then solve Eqs. (26) and (31) simul- 
taneously for 6(0a/dc). 
mm a a Lin and the author have calculated the velocity along 
~ the axis of symmetry behind a detached shock wave as 
a power series expansion from the nose of the body.! 
at! For a body of revolution whose curvature at its nose is 
ko, the expression becomes 
Fic. 1. Values of g and h. 
= ax + bx? +... (32) 
4 + In the above, a represents an undetermined constant, 
| ~ and x represents the distance from the nose of the body 
A i. to the point in question. In this coordinate system, the 
shock is located at = —6 and the nose of the body at 
A | | x = 0; bis related to a through, 
“2 ~2.0 b = 2aky — \F(a/da)?* (33) 
| where 
| 
= 4y(1 — — (vy — 1)] X 
2+ (y¥- 4) 
4 -1.0 1 um y+ - 
2 3 4 : ( M. ‘) = 
M, 


Fic. 2. Values of A and yz. 


= g/u (29) 


This, in effect, is approximating the true velocity dis- 
tribution by a straight line that starts at the known 
value of wu: and crops to zero at the nose of the body. 
The slopes of tle approximate and true velocity dis- 
tributions are equal at the shock. 

If one uses three terms of the series, one will have a 
quadratic that can be solved for 5(0a/0e) to give 

= (u + Wu? — 2gh)/h (30) 

The povitive sign is used in front of the radical. 

This in effect is approximating the true velocity dis- 
tribution by a parabola that starts at the known value of 
uz and drops off to zero at the nose of the body. The 
(Ou/Ox)2’s and the (0?u/0x?)2’s of the true curve and the 
approximation are equal at the shock. 

If one includes the fourth term of the series, one will 
have two unknowns—namely, 6(0a/0c) and 6°(0'a/ 
0o*)—-since (03u/dx*). is known to consist of a (0a/dc)* 
and a 0’a/de* term. To solve for 6(0a/0c), therefore, 
another equation is required. This is obtained by dif- 
ferentiating Eq. (24), applying it at the nose of the 
body, and multiplying through by p2%25/po. This re- 


(35) 
pole Po 


Values of \ for a Mach Number range of 1.25 to 4.50 
are shown plotted in Fig. 2. 


By differentiating Eq. (32) once and twice, respec- 
tively, and then applying conditions at the nose of the 
body (x = 0 in this case), one obtains expressions for 
(du/dx)y and (d?u/dx*)y in terms of a and b. It should 
be noted here that the inclusion of higher order terms in 
Eq. (32) will not affect these two derivatives of velocity 
or the a — 6b relationship, Eq. (33).. Therefore, any 
inaccuracy in this development will result from insuffi- 
ciently high order terms in the expansion from the shock, 
Eq. (24), and not from the expansion from the nose of 
the body, Eq. (32). 

If Eq. (24) is differentiated twice, applied at the nose 


of the body, (x = 6), and multiplied by p2%26°/po, one 
obtains 


= ki —{—]}63+... (36) 
po \dx?/o Oa Ox*/ 
One now proceeds to eliminate the (du/dx)o and 
(d?u/dx*)y in Eqs. (31) and (36) by replacing them by 
their values in terms of a and b. By the use of Eq. 
(33), both Eqs. (31) and (36) can be written in terms of 
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the three variables, 6(0a/0e), a, and (0°u/0x*)26°. 


a new expression that. involves only 6(0a/0c) and 


The 
ais then eliminated between the two equations to give 
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(0%u/Ox*).6*. This new expression can now be solved 
simultaneously with Eq. (26) to give a quadratic equa- 
tion whose roots are 


Oa — 3) = Vyu2(4kod — 3)? — 12¢(2kod — 1) [h(kod — 1) + AJ 


(37) 


Oo —2[h(kod — 1) + d] 


The positive sign is used in front of the radical. 

This, in effect, is approximating the true velocity 
distribution by a third degree curve that starts at the 
known value of uw: and drops to zero at the nose of the 
body. The (Ou/Ox)2’s and (0?u/0x)2’s of the true curve 
and the approximation are equal at the shock, and an 
exact relation, which is satisfied by the true (du/dx)o 
and (d*u/dx*), at the body nose, is also satisfied in the 
approximation. 

Asa rough check on the theoretical values of 6(0x/0c), 
some experimental values of 5(0a/0c) were obtained 
from photographs of spheres in high-speed flight.* 
In these experimental values, the distance 6 was ob- 
tained by measuring the distance from the outermost 
dark shadow of the shock to the rearmost point on the 
sphere and then subtracting from this distance the 
diameter of the sphere obtained from the uppermost 
and lowermost points on the sphere. This assumes that 
all the distortion in the photograph of the sphere caused 
through refraction by the shock occurs at the front of 
the sphere. The O0a/0o of the shock was determined 
by measuring the radius of the circular are that best 
fitted the points of the shock near the axis of sym- 
metry. One should bear in mind that these experi- 
mental values are crude because of the distortion. and 
the rough method of obtaining the data. 

A more accurate determination of the shock is fur- 
nished by Ladenburg and others? by the interferometric 
method. They observed the steady supersonic flow 
past a sphere at a Mach Number of 1.7. Their value 
of 6(0a/Oc) is also reproduced in Fig. 3. It is seen that 
the agreement is excellent. 

The values of 6(0a/0c) obtained through the use of 
Eq. (29) and (30) are plotted in Fig. 3, together with 
the experimental data. The theoretical computations 
seem to be slightly higher than the trend of the rough 
experimental data. Also, the parabolic approxima- 
tion, Eq. (30), seems to be further off than the linear 
approximation, Eq. (29). This can be explained, how- 
ever, by the following brief qualitative analysis of the 
nature of the velocity distribution curve. 

From an investigation of the expansion from the nose 
of the body, Eq. (32), it is evident that the coefficient a 
must be negative, since for small negative distances of 
* the velocity must be positive quantity. The sign of 
b cannot, however, be determined directly from Eq. 


* Photographs by Dr. A. C. Charters, loaned to the author by 
Prof. H. W. Emmons, of Harvard University. 


(33), since the 2ak) term is negative whereas the 
—\F(0a/dc)? term is positive. But an indication of 
the sign can be obtained by assuming the axially sym- 
metric value of } to be of the same sign as that in the 
two-dimensional case. The power series gave, in the 
two-dimensional case, b = */2ako. Since ko is positive 
and a was found to be negative, } will also be negative. 
With these signs of a and b, (du/dx)o and (d*u/dx*)o 
at the body nose are both negative. To investigate 
conditions at the shock, Eq. (26), with the 6’s replaced 
by x’s, isexamined. In this equation, x = 0 represents 
conditions at the shock. It is evident that, since yu 
and 0a/dc are both negative, the (du/dx)2 will also 
be negative. Also, since g and h are both positive, the 
uz and (d?u/dx*), will be positive. If one gathers to- 
gether this information on the conditions of velocity at 
the body nose and at the shock and then draws the sim- 
plest possible curve, one obtains an S-shaped curve— 
i.e., one with a point of inflection between the body 
nose and the shock. Fig. 4 illustrates this point and 
also how the linear and parabolic approximations fit 
such a curve. 

The cubic approximation, Eq. (37), was used to 
check the five experimental values. In these checks, 
one must remember that the measured values of ko and 
6 for use in Eq. (37) were crude. A table of the theo- 
retical and experimental values is presented in Fig. 5. 
The theory seems to have broken down for M, = 1.53. 
This can be explained by investigating the radical of 
Eq. (37). Because the second term inside the radical 
is multiplied by the quantity (2k96 — 1), any kod that 
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~TRUE CURVE 
7—PARABOLIC APPROXIMATION 


LINEAR APPROXIMATION 


uy 


0 
SHOCK NOSE 
x 


Fic. 4. Velocity distribution behind the shock. (Not to scale). 
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Fic. 5. Agreement between experiment (Charter’s Results) and 
the theoretical cubic approximation | Eq. (37)]. 


approaches 0.500 will make the second term extremely 
small, thereby making the whole numerator approach 
zero. Physically, this means that, for a 6 approaching 
a half of the radius of curvature of the body nose, the 
interval over which the power series expansion from the 
shock is taken becomes too large and the whole method 
herein described breaks down. For the other points 
where the value of ko6 was not close to 0.500, the theo- 
retical results seem to have agreed favorably with the 
crudely determined experimental data. 

Eq. (37) then seems to give a good estimate of the 
6(0a/0c), introducing the curvature of the body nose 
into the problem. Its accuracy can be improved of 
course by introducing higher order terms into the series, 
Eq. (24). 


(4) Serres EXPANSION AROUND THE NOSE OF THE 
Bopy 


If the velocity distribution is approximated by a 
power series from the nose of the body rather than from 
the shock, it becomes possible to obtain 6 and 0a/d¢ 


individually, given only the curvature of the body nose, 
ko, and the free stream Mach Number, M4. 


As was stated in the previous sectién, the velocity 
behind a shock can be expanded as a power series from 
the nose of the body, Eq. (32). If this expansion js 
applied to conditions at the nose of the shock (x = 
—6), an equation involving the three unknowns 
0a/do, a, and 6 results. Two more relations involving 
these same unknowns are obtained by differentiating 
Eq. (32) once and twice, respectively, and again setting 
in conditions at the shock. The three equations ob- 
tained are: 


1 


V bo, U2 gF ad + 2akod 


V po/ po dx 2 Oa 


Laks + arr (52) ‘+... 
Oo 
V bo/ po dx? 2 Oa 


Eqs. (38) can be brought into a simpler nondimen- 
sional form by making the following substitutions: 
n = a/ko 
é = kod 
j= 


This changes Eqs. (38) to 


gF = + 2né? — +... 
—pF(SE) = — + +...7 (40) 
hF(¢é)? = — 


By eliminating first the né? term and then the 7, one 
obtains 


(38) 


(39) 


De 


The above is identical to Eq. (30), the parabolic 
velocity approximation matching the (du/dx), and 
(d?u/dx*)2 at the shock. 

Continuing the solution of the three equations results 
in 


= @a/de) = (43) 


The above equations completely solve for the two 
parameters 6 and Oa/0c. Unfortunately, the results 
obtained from Eqs. (42) and (43) disagree greatly with 
the experimental data. For an M, = 1.95, Eqs. (42) 
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and (43) yield € = 0.007 and ¢ = —35.0. This con- 
trasts markedly with the experimental values of & = 


0.37 and ¢ = —0.44. This disagreement was more or 
less expected, since the power series expansion from the 
body nose which goes up to x? only is a crude approxi- 
mation with which to match a curve that should prob- 
ably have a point of inflection over its range. If Eq. 
(32) were developed to include an x* term and then ap- 
plied to the conditions at the shock, much better agree- 
ment with experiment is presumed. 


(5) CONCLUSIONS 


It seems that the power series expansions are actually 
useful for such calculation in the detached shock wave 
problem and that further work should be done in this 
direction. The prime necessity now seems to be the 
calculation of higher order approximations for more 
accurate results. This means, basically, the evaluation 
of (0'u/Ox*)2 and (0‘u/Ox*)2 and also the expansion of 
the power series to include x* and x‘ terms. 

An interesting result that this investigation seems to 
bring out is the probable shape of the velocity distribu- 
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tion curve between the shock and the body. Because 
of this S-shaped curve, it appears that both the linear 
velocity approximation, Eq. (29), and the parabolic 
velocity approximation, Eq. (30), which are both inde- 
pendent of the body nose curvature, are too crude to 
give anything beyond a rough estimation of the quan- 
tity 6(0a/0c). The cubic velocity approximation, 
Eq. (37), is more refined. It should, however, be 
checked against more accurate experimental data to 
determine whether or not it is capable of giving useful 
engineering results. 
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Errata 


Becker, H.,* The Optimum Proportions of a Long Unstiffened Circular Cylinder in Pure Bending 
Journal of the Aeronautical Sciences, Vol. 15, No. 10, pp. 616-624, October, 1948. 


(1) The third equation after Eq. (9), page 619, should read 
= — YD 
instead of 
= — 
(2) The expression following Eq. (10), page 619, should be revised from 
R/Ro = ..... = (00/0)/(Wo/W) 
to 
Ri = (00/0) /(W/Ws) 


(3) In Fig. 5, page 620, the table of o:/o2 vs. n should be altered as shown: 


Present Value Corrected Value 
of n of 
1.048 12 1.084 
1.084 20 1.048 


(4) The comparison of weight errors for several materials, using the various plastic modulus theories, is mis- 
leading in its present form (Table 1, page 622, of the original article). The tabulated values of R/Ro (and, sub- 
sequently, of W/Wo) were obtained by entering Fig. 5 with ¢/a9 = oyp/o0. This implies the use of o,, in the error 
equations [Eqs. (10) through (11)], together with the corresponding effective moduli. However, the use of 
Cyp for o; in the simplified design equation [Eq. (16), page 623] disregards the question of effective moduli. Conse- 
quently, the weight and radius accuracies of the simplified design equation exceed those of the error equations, 
wherein the radical changes in E/E, affect the errors appreciably. 

In order to test the accuracy of Eq. (16), the following development may be used: 


Ro = Kr(ME/oy,*)'/* (16) 
or 
Reyp/Rer = (El) 
using Reyp = Ro and Rey = Ry in Eq. (16), and 
Reyp/Ro = (Reyp/Rei)(Rei/Ro) (E2) 


The values of Rei/Ro may be found, from Fig. 3, to be close to unity for all o;/o2, as is shown by the weight 
error column. Furthermore, the largest weight error, using Reyp/Ro = Reyp/Roi, amounts to only 1.5 per cent 
(for type 304 stainless-steel sheet, for which o; = 30,600 lbs. per sq.in., ¢yp = 36,700 Ibs. per sq.in., Reyp/Rei = 
0.834 using Eq. (E1), and (W/Wo) — 1 = 1.5 per cent using Fig. 5), while the weight errors for most of the ma- 
terials are negligible. This means that the simplified design procedure using the tangent modulus theory will 
deviate only 2 to 3!/2 per cent from optimum. The large weight discrepancies obtained by using cy», together 
with the tangent modulus and reduced modulus theories in the error equations, imply the inadequacy of any but 
the secant modulus theory for simplified design. This is disproved by the preceding analysis, which indicates 
that, no matter which plastic bending theory is correct, the simplified design equation (with the appropriate 
value of Kr) may be employed with little concern for weight error. Its use is not restricted to the secant modulus 
theory as was indicated in the article. 


* EDO Corporation, College Point, N.Y. 
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Bending and Buckling of Sandwich Beams’ 


N. J. HOFF? ann S. E. MAUTNER? 
Polytechnic Institute of Brooklyn and Skydyne Inc. 


SUMMARY 


Formulas are derived for the calculation of the deflections of 
sandwich beams and of the buckling load of sandwich columns. 
Beam and column tests carried out on 64 specimens are de- 
scribed. The agreement between theory and experiment was 
found to be good in the column tests when the core was balsa- 
wood and in all the beam tests. A number of sandwich columns 
having a cellular cellulose acetate core failed under loads smaller 
than those predicted by theory. The early buckling was 
probably caused by weak spots in the core material or the 
adhesive. 


INTRODUCTION 


_— DEVELOPMENT OF sandwich-type structural 
elements for airplane construction has recently 
received considerable attention. This is because the 
sandwich skin remains smooth under the heavy loads 
to which high-speed aircraft is subjected. It is known 
that the smoothness of the surface of the wing and the 
fuselage and the true shape of the contours must be 
maintained more and more accurately as the high 
speed of the plane is increased. However, the thin- 
walled aluminum-alloy structures used almost exclu- 
sively in airplane construction at the outset of World 
War II did not satisfy this requirement. Their thin 
skin buckled or wrinkled under comparatively light 
loads. 

To eliminate this type of instability, imaginative 
designers introduced the sandwich skin consisting of 
two very thin layers of high-strength ‘‘face’’ material, 
between which a thick layer of ultra-lightweight “‘core”’ 
is sandwiched. Materials available at present for the 
faces are: plywood, papreg, aluminum alloys, stainless 
steel, and fiberglas; and for the core: balsawood, 
expanded cellulose acetate, polystyrene, and rubber, 
as well as built-up grids of plywood, synthetic ma- 
terials, and metal. The moment of inertia of the cross 
section of the sandwich skin is large because of the 
comparatively great distance between the two faces, 
and thus the buckling formulas derived from classical 
plate theory give buckling stresses that, as a rule, 
far exceed the yield stress of the face material. 

The classical formulas, however, are not directly 
applicable to the calculation of the instability of the 
sandwich skin, since the lightweight core permits 
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unusually large shearing deformations. Moreover, 
relative displacements of the two faces are possible 
because of the small extensional rigidity of some of the 
core materials used. With the aid of the principle of 
virtual displacements, a new theory of bending and 
buckling of sandwich-type structural elements is 
developed in the present paper from a consideration 
of the essential portions of the strain energy stored 
during deformations. These are, in the opinion of the 
writers, the strain energy caused in the faces by exten- 
sion and bending and that caused in the core by shear 
and by extension perpendicular to the plane of the 
faces. Quantities disregarded are, therefore, the strain 
energy of shear in the faces and that of extension 
parallel to the faces in the core. 

The analysis yielded comparatively simple formulas 
that can be used in routine calculations. Beam and 
column tests carried out on 64 specimens (see ‘Test 
Specimens” tables) proved the general validity of 
the theory. However, the failing loads of a number 
of sandwich columns were considerably below the 
values predicted when the core material was cellular 
cellulose acetate. It is recommended, therefore, that 
conservative estimates be made of the modulus values 
whenever expanded materials are used for the core. 


RESULTS OF THE ANALYSIS 


The analysis was carried out for a sandwich beam 
whose core is weak as compared to its faces so that the 
bending moment is resisted by the faces alone. This 
assumption is satisfactory when the core is an expanded 
high polymer, a grid (honeycomb construction), or 
balsawood with its grain perpendicular to the faces. 
When the grain of a wood core is arranged parallel to 
the axis of the beam, the contribution of the core to 
the moment equilibrium cannot, as a rule, be neglected. 


Deflections 


When the cantilever shown in Fig. 1 is subjected to 
a concentrated lateral load W, the average end de- 
flection of the two faces is 


d, = (¥/3)(WL*)/(El), (1) 
where 
= 1+ (la) 
and 
C = [3/(pL)*)[1 — (tanh pL/pL)) (1b) 


Moreover 
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= (Gb)/(EI)* (2) 


(EI)* = (ED) + (2ET),| (3) 


and the moment of inertia quantities are defined as 


(EI); = (1/2)t(6 + t)?E, (4a) 
(2ED), = (1/6)tEy (4b) 
(ED, = (ED: + (ED, (4c) 


where /, is Young's modulus for the face; G, the shear 
modulus for the core; W, the load; L, the length of 
the cantilever beam; 1¢, the thickness of one face; 3, 
the thickness of the core; and the width of the beam 
is taken as | in. 

Values of C are plotted against pl in Fig. 2. The 
following particular cases deserve attention: 

(a) When pL < 0.1, 


dy = (WL')/3(2ED, (5) 


Consequently, when the shearing rigidity of the core 
is extremely small, the sandwich beam deflects like a 
beam consisting only of the two faces. 

(b) When pL > 100, 


(6) 


Since in most cases (//), differs but little from (£/),, 
the deflection is simply the sum of the bending de- 
flection and the shearing deflection when pL > 100. 

(c) When pL (because G— while ¢ remains 
finite), 


d, = (WL*)/(3EI,) (7) 


The deflection is the simple bending deflection of the 
Euler-Bernoulli theory when the core is infinitely rigid 
in shear. 
The mid-deflection of a beam on two simple supports 
is 
d, = (y/48)(WL?)/(ED, (8) 


The symbols have the same meaning as before. 

Simultaneously with the deflection of the free end 
of the cantilever, the two faces of the sandwich beam 
approach each other because the weak core is com- 
pressed. Although this decrease d; in the distance 
between the faces at the free end of the cantilever is 
insignificant in sandwich beams of normal design, the 
formula is given here: 


ds = (¢/x)(Wb)/(E.L) (9) 


where /, is Young’s modulus of the core for compres- 
sion perpendicular to the faces, and 
DL(sinh 2DL — sin 2DL) 


= % 
mi cosh? DL + cos? DL (9a) 


where D can be calculated from the equation 
D* = (3x/bt*)(E./Ey) (10) 


In most calculations x can be assumed to be 2. More 
is said about x under ‘Calculation of the Strain 


Energy.” 
Whenever DL > 3, Eq. (9) can be simplified to read 
ds = (WbD)/(E.x) (11) 

Buckling 


When a sandwich column is pinned at its two ends 
and has a core whose modulus, in the direction of the 
axis of the column, is small enough to justify the as- 
sumption that the entire compressive load is carried 
by the two faces alone, then the buckling load can be 
calculated from the formula 


Per = Pi{ + + + Gb) (12) 


where 
P, = w(EI),/L? (12a) 
P, = 2°(2EI),/L? (12b) 
and the bending rigidity quantities were defined in 
Eqs. (4). 


The following particular cases are of interest: 
(a) When the bending rigidity (2£/), of the faces 
is negligibly small, Eq. (12) reduces to 


P., = P,Gb/(P; + Gb) (13) 


As may be seen from Timoshenko’s Theory of Elastic 
Stability,! this buckling formula is essentially the 
same as the one obtained by Engesser when he investi- 
gated the effect of the shearing deformations upon the 
buckling load of ordinary columns. 
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(b) When the shearing rigidity Gb is large, Eq. 
(12) becomes 


P., = + (14) 


This is the ordinary Euler formula for the sandwich 
column, neglecting the capacity of the core to carry 
compression. 

(c) When the shearing rigidity Gb is negligibly 
small the buckling load is 


= (15) 


This is the Euler formula for the two faces that are 
independent when the shearing -igidity of the core 
vanishes. 

When the ends of the column are fixed, P; and P; 
in Eqs. (12a) and (12b) must be multiplied by 4. 
When the buckling stress obtained from Eq. (12) 
exceeds the proportional limit, the value of E in Eqs. 
(12a) and (12b) must be replaced by £,, the tangent 
modulus value corresponding to the buckling stress 
obtained from Eq. (12). To find the appropriate 
value of the modulus beyond the elastic limit, a trial- 
and-error process is recommended. 

The faces can also buckle in the form of a short 
wave ripple when the shearing and extensional rigidities 
of the core are low. The buckling load expressions 
obtained for this type of instability are somewhat 
involved. Since, however, they are fully equivalent 
to those given in reference 2, the simplified formula 
recommended there is quoted: 


Ser (1/2) VE,E.G, (16) 


where f,, is the critical compressive stress acting on the 
faces, E, is Young’s modulus of the faces parallel to 
the axis of the beam, £, is Young’s modulus of the 
core perpendicular to the faces, and G, is the shear 
modulus of the core associated with the directions of 
the axis of the beam and of the perpendicular to the 
faces. 

This wrinkling or ripple type of instability is inde- 
pendent of the end conditions of the beam. It is 
substantially a buckling of the face layer restrained 
by the action of the elastic foundation provided by the 
core. Buckling according to Eq. (16) will occur if the 
corresponding buckling load is smaller than that 
calculated from Eq. (12), and vice versa. 


EXPERIMENTS AND COMPARISON OF THEIR RESULTS 
WITH THEORY 


Compression Tests of Cellular Cellulose Acetate Cubes 


Five cubes of 3-in. edge length, composed of six 
layers of cellular cellulose acetate glued together, 
were tested in compression on a Fairbanks lever-type 
testing machine of 2,000 Ibs. capacity. (See Fig. 3.) 
The machine was only used to measure the load, and 
the uniform compressive strain was realized by dis- 


Fic. 3. Compression test. 
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Fic. 4. Compression test of CCA cube specimen D. 


placing a steel plate resting upon the cube by turning 
four adjusting screws. The displacements were meas- 
ured by four Ames dial gages having 0.001 or 0.0001 
in. subdivisions. Average displacements were plotted 
against loads. Fig. 4 shows the plot obtained with 
specimen D. Young’s modulus for the material was 
calculated from the slope of the experimental line. 
It varied from about 4,000 to about 5,000 Ibs. per 
sq.in. The material weighed 4.9 Ibs. per cu.ft. 
Changes in the length of the specimen in the trans- 
verse direction were measured by four more Ames 
dial gages. In most cases slight elongations were 
observed, but in a few others a shortening occurred. 
The measured values were plotted against the load, 
and the average value of Poisson's ratio was determined 
as the ratio of the average slope of the transverse elon- 
gation curves to the slope of the plot of the shortening 
in the direction of the compressive load. Values ob- 
tained for Poisson's ratio varied between 0.04 and 0.14, 
with an average of 0.09. A typical example of the 
transverse deformation curves is shown in Fig. 5 
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Fic. 6. Bending test. 


On the basis of these tests it appears permissible to 
assume Poisson’s ratio as zero in the calculations. 
Consequently, the shear modulus may be taken as one- 
half of Young’s modulus. 


Bending Tests of Sandwich Beams 


The setup for the bending tests may be seen in 
Fig. 6. The simply supported beam was loaded in 


CORE THICKNESS, IN. 


Fic. 7. Comparison of calculated and measured beam deflec- 
tions. Length between supports, 10 in. 


the middle by dead weights, and the deflection of the 
midpoint was measured by means of an Ames dial gage. 
The purpose of these tests was twofold: first, to 
determine the shearing rigidity of the core material 
with the aid of the theoretical deflection formulas; 
and, second, to check the theory by comparing its 
predictions with the deflections observed in tests with 
beams of various lengths and various }/t ratios. The 
width of the beams was | in. 

Each beam was loaded several times, and the de- 
flections measured were plotted against the loads. 
The result was always a straight line, the slope of which 
was computed. In Fig. 7 the slopes obtained with 
six specimens of 10-in. supported length and having a 
transverse balsa core (the fibers of the balsa core are 
perpendicular to the planes of the faces) are plotted 
against the thickness 6 of the core. The figure also 
contains results obtained in similar tests carried out 
with seven specimens having a cellular cellulose acetate 
core that weighed 4.9 lbs. per cuft. In Figs. 8-10 
thé ordinate is the slope of the deflection curve, while 
the abscissa is the distance between supports. The 
specimens corresponding to the first of these figures 
had a '/2-in. thick transverse balsa core weighing 5 
Ibs. per cu.ft.; those corresponding to the second and 
third figure had a '/,-in. thick cellular cellulose acetate 
core weighing 5.5 Ibs. per cu.ft. The faces were of 
0.032-in. thick 24ST Alclad in Fig. 9 and of 0.081-in. 
thick 52 S-1/2H aluminum alloy in Fig. 10. 
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Fic. 8. Comparison of calculated and measured beam deflec- 
tions. 


In Figs. 7 and 8 good agreement between experi- 
ment and theory was found when the shear modulus 
of the balsa wood was taken as 24,000 lbs. per sq.in. 
The shear modulus of the cellular cellulose acetate 
(Fig. 7) was assumed as 2,500 Ibs. per sq.in., which is 
about one-half the value of the Young’s modulus found 
in the compression tests. 

The cellular cellulose acetate of the specimens repre- 
sented in Figs. 9 and 10 was heavier and of higher 
mechanical properties than that used in the earlier 
tests. Theoretical curves consistent with the experi- 
mental points were obtained when the shear modulus 
was assumed as 5,000 Ibs. per sq.in. 

Deflections were calculated from Eqs. (1) and (9). 
The contribution of the latter—that is, the flattening 
of the core—was always insignificant. The values of 
the parameter PL varied between 2 and 75. This is 
the range in which the interaction between shearing 
and face bending is pronounced, and thus the simplified 
formulas are not applicable. 


Buckling Tests of Sandwich Columns 


Most of the specimens first tested as beams below 
the elastic limit, as well as an additional number of 
specimens of similar construction, were tested to 
failure in axial compression. They were placed with 
flat ends in a lever-type Tinius-Olsen testing machine 
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Fic. 10. Comparison of calculated and measured beam deflec- 
tions. CCA core; 0.081-in. face. 


of 2,000—20,000-Ib. capacity, and the lateral deflections 
of the midpoint were measured by means of an Ames 
dial gage. (See Fig. 11.) The purpose of the lateral 
deflection readings was to obtain a check on the uni- 
formity of loading and the straightness of the specimen. 
Experiments carried out by the Aluminum Company 
of America with round rods and structural shapes*® 
and earlier tests at PIBAL, as well as the shape of 
deflections observed in the present test series, indicate 
that the “flat end’’ condition is practically equivalent 


| 
| 

he 

ze. 

to 

ial 

1a y, DD 

its 

th 

he 

le- 

ls. 

ch 

th 

a 

= 

od 

it 

te 

0 

le 

| 

5 

d 

e 

yf 

1. 


712 JOURNAL OF THE AERONAUTICAL SCIENCES—DECEMBER, 1948 


= 


Fic. 11. Column test. 
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to rigid end fixation in the case of the present test 
specimens. 

In Figs. 12 and 13 are plotted the buckling stresses 
of the specimens having a balsa wood core. In Figs. 
14 to 17 are plotted those of the specimens having a 
cellular c-llulose acetate core. In Fig. 12 the abscissa 
is the thickness of the core, and in the other figures the 
abscissa is the length of the column. The faces of the 
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specimens of Figs. 12 and 13 were of 0.040-in. thick 
24ST Alclad, the faces of the specimens presented in 
Figs. 14 and 15 were of 0.032-in. thick 24ST Alclad, 
while the specimens of Figs. 16 and 17 were constructed 
with 0.081-in. thick 52 S-1/2H aluminum alloy. 

Buckling always occurred according to a full sine 
wave pattern corresponding to rigid end fixation, and 
a considerable portion of, the deflections was visibly of 
the shear type. With most specimens the deflection 
of the midpoint of the column was small until about 
one-half the buckling load was reached. When larger 
than average deflections were observed in the early 
stages of loading, the specimen usually failed under a 
small load. 

At failure, a transverse crack developed all across 
the core, usually near a quarterpoint of the length of 
the column, and either above or below the crack the 
core separated from one face over the entire distance 
extending from the crack to the end of the column. 

The only column having a 1-in. thick core buckled 
sideways—that is, the displacements occurred in 
planes parallel to the faces, not perpendicular to them 
as in all the other specimens. 

The theoretical buckling load was calculated from 
Eq. (12). The shear moduli used in the computations 
were those determined earlier from the bending tests— 
that is, G = 24,000 Ibs. per sq.in. for the balsa core of 
the specimens in Figs. 12 and 13 and G = 5,000 Ibs. 
per sq.in for the celiular cellulose acetate of the speci- 
mens in Figs. 14 to 17. The buckling stress of the 
specimens having a balsa core was higher than the 
elastic limit. For this reason the stress values were 
corrected using the Engesser-von K4rmaén_ theory,‘ 
as well as the tangent modulus theory. The formula 
used for the reduced modulus was 


Ena. = 2EE,/(E + E,) (17) 


The values of the 


where F; is the tangent modulus. 
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tangent modulus were taken from a publication of the 
A'uminum Company of America.* 


The agreement between theory and experiment is 
satisfactory with specimens having a balsa core. 
When the core is cellular cellulose acetate, the scatter 
of the experimental points is considerable. The 
addition of the second theoretical curve—namely, the 
one corresponding to a shear modulus of 2,500 Ibs. 
per sq.in.—may help to explain this scatter. The 
effect of a change in the shear modulus upon the 
buckling stress is small when the slenderness ratio of 
the sandwich column is either extremely large or 
extremely small. In the former case the buckling 
stress is practically equal to the buckling stress obtain- 
able from the Euler formula, using the moment of 
inertia of the column, and in the latter to that cal- 
culated from the Euler formula in which the moment 
of inertia of the faces is substituted. In the inter- 
mediate range, however, a reduction in the shear 
modulus value results in substantially diminished 
buckling stresses. 


Expanded core materials are never complete'y uni- 
form. In particular, when they are cut down to small 
thicknesses, the random distribution of the pores can 
easily result in considerable variations in the rigidity 
of different cross sections. When the shear modulus 
of the core is determined from a beam test, every 
cross section is of equal importance and the, palin 
value obtained is an arithmetic mean. Information, 
test, however, regions of large shearing de ‘- smalle« 
can shift to regions of small shearing rigidi’ 
the apparent shearing rigidity of the cote is.) 
in a column test than in a beam test, and the difference 
between the two values depends upon the random 
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variations of the rigidity from one cross section to the 
other. 


On the basis of these considerations the agreement 
between theory and experiment appears to be reason- 
able in Figs. 14 to 17. The average experimental 
values, including the balsa core specimens, are about 
80 per cent of those predicted by theory. A few 
individual specimens, poorly constructed or centered, 
failed at about one-half the theoretical buckling load. 
It seems advisable, however, to estimate the shearing 


.tigidity of expanded cores conservatively when the 


buckling load of a sandwich column is predicted with 
the aid of the theory. 


DEVELOPMENT OF THE THEORY 


Calculation of the Strain Energy in a Sandwich Beam 


A cantilever sandwich-type beam of unit width is 
shown in Fig. 18. If an end load is acting upon it, 
it is assumed that the load 2ot is applied to the face 
material only, since the load-carrying capacity of the 
core is small as compared to that of the faces. As 
long as a stability limit is not exceeded, the end load 
will only cause a uniform axial compression. This 
compressed state of equilibrium is considered as the 


initial state of the beam when the variation of the 


is undertaken. 
he defle 
th ey represeted shape is defined by the three functions 


47, “are ot the coordinate x. As shown in Fig. 18 


extensional and shearing deformations 


anu compression of the core, respectively. The strain 
energy stored in the beam can then be calculated as 


follows: 
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The strain energy of extension U,,,. stored in the two 
faces is completely defined by the axial displacement 
function u: 


L 
= edV = u’? dx (18) 


where F, is the modulus of elasticity of the face ma- 
terial. 

The strain energy of bending U, depends upon both 
v and w, and the two are combined differently for the 


upper and lower face: 


= (1/2)(El), (v” + w")*dx + 


L 
(1/2)(ED), (v" — w")*dx (19) 
0 


where 


(EI), = (1/12)E,#° (19a) 


is the bending rigidity of one face alone. 

If the displacement u is assumed to prevail at the 
middle of the thickness ¢ of the face, then the angle of 
shear y; caused by it in the core is 


= 2u/(b + (20a) 
The angle of shear y2 due to 7 is 
= —v’ (20b) 


Since « and v represent antisymmetric and w repre- 
sents symmetric distortions, the strain energy of shear 
is composed of two portions, one due to the antisym- 
metric distortions and the other due to the symmetric 


ones. The former can be written as 


(v1 + ¥2)"dx = 
0 


L 
(1/2)Gb {[2u/(b + — v’}%dx (20) 
0 : 


where G is the shear modulus of the core. 

The shear strain energy caused by w will not be 
considered separately but only in conjunction with the 
strain energy of compression of the core. In_ this 
calculation it is assumed that there is a proportionality 
between the compressive loads p applied to the two 
faces at any value x of the axial coordinate and the 
displacement w at x. In reference 2 such a propor- 
tionality was proved to exist in the case when the load 
p is distributed according to a sine law: 


p = posin (rx/L’) (21a) 


The deflectian can be calculated from an equation 
developed in reference 2 and rewritten here with slight 


changes in notation: 
b = xE(w/b) 


where E, is the modulus of elasticity of the core and 


(21b) 
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4(rb/2L’) cosh? (rb/2L’) 


21c) 
X= 3 cosh (wb/ 2L’) sinh (xb/2L') — (xb/2L') 


Values of x are plotted in Fig. 19 against the ratio of 
wave length to core thickness. It can be seen from the 
figure that x = 2 = constant in good approximation 
when the wave length L’ is greater than the core thick- 
ness b. The assumption of proportionality is justified, 
therefore, in the following two cases: (1) when the 
distortions are sinusoidal; (2) when the distortions are 
nonsinusoidal and the coefficients of a Fourier expan- 
sion of the distortions are negligibly small for all those 
terms that correspond to wave lengths smaller than the 
core thickness. 

It is noted that in the derivation of Eqs. (21b) and 
(2lc) Poisson’s ratio was assumed to,be zero for the 
core. 

If the associated shear strain energy is included, 
the strain energy of compression UL’, stored in the core 
is equal to the work done by the distributed loads p 
during the deflections of both faces: 


L 
U,= pw dx = x(E,/b) dx = (21) 
0 0 


The total strain energy UL’ is the sum of the four 
quantities given: 


U = Vert, + U, T Us + U, (22) 


The Potential of the Loads 
The potential V, of the distributed load can be 


written as 
L 
V, = -f q(v + w)dx (23) 
0 
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Fic. 19. Compression w of core under unit load p. Sinusoidal 
distribution with wave length L’. 


The. potential V, of the concentrated load is 
Ve = + (24) 


The potential V, of the end load 2ot is the negative 
product of the load and the shortening of the distance 
L because of bending: 


V, = —(1/2)ot X 


L L 
Lf (v’ + w’)*dx + f (v’ — | (25) 
0 


The potential of all the external forces is the sum of 
the quantities expressed in Eqs. (23)—(25): 


V= Vat Vet V, (26) 


Equations of Equilibrium of the Cantilever when q = « = 0 


The equations of equilibrium can be derived from the requirement that the first variation of the total potential 


must vanish (principle of virtual displacements) : 


+ V) =0 


(27) 


L 
V) = + (El) + w") + dw”) + (ED — w")(6v" — bw”) + 
0 


Gb{[2u/(b + t)] — { [26u/(b + — + 2x — + = 0 (28) 


With 


when x = 0, and 


when x = L, integration by parts yields 


(29a) 


(29b) 
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TEST SPECIMENS 
Core Length of Buckling 
Thickness, Specimen, Tested in ’ 
Specimen Core Material In. Face In. Bending Lbs. 
2-A-1 Balsa 1 24ST Alclad 12 rae 3,510 
2-A-2 (Transverse) 1 0.040 in. thick 12 = ee 
2-A-3 5 lbs. per cu.ft. 1 0.040 in. thick 12 x dt 
2-B-1 5 Ibs. per cu.ft. 1/, 0.040 in. thick 12 x pee 
2-B-2 5 lbs. per cu.ft. 1/, 0.040 in. thick 12 x rae 
2-B-3 5 Ibs. per cu.ft. 1/s 0.040 in. thick 12 — 3,110 
2-C-1 5 Ibs. per cu.ft. V/, 0.040 in. thick 12 ee 1,840 
2-C-2 5 Ibs. per cu.ft. 1/, 0.040 in. thick 12 x ide 
2-C-3 5 Ibs. per cu.ft. 1/, 0.040 in. thick 12 x —- 
3-A-2 5 lbs. per cu.ft. 1/s 0.040 in. thick 10 x 3,270 
3-B-1 5 Ibs. per cu.ft. 1/, 0.040 in. thick 15 x 2,840 
3-C-1 5 Ibs. per cu.ft. 1/, 0.040 in. thick 20 x 2,190 
3-D-] 5 Ibs. per cu.ft. 1/, 0.040 in. thick 25 x 2,120 
2-D-2 CCA 1 0.040 in. thick 12 x bee 
2-D-3 (Cellular cellulose acetate) 1 - 0.040 in. thick 12 x sed 
2-E-1 4.9 Ibs. per cu.ft. 1/9 0.040 in. thick 12 x ea 
2-E-3 4.9 lbs. per cu.ft. 1/9 0.040 in. thick 12 x = 
2-F-1 4.9 Ibs. per cu.ft. 1/, 0.040 in. thick 12 x we 
2-F-2 4.9 lbs. per cu.ft. 1/, 0.040 in. thick 12 x ae 
2-F-3 4.9 lbs. per cu.ft. 1/, 0.040 in. thick 12 x os 
4-B-1 5.5 lbs. per cu.ft. 1/, 0.032 in. thick 10 x 525 
4-B-2 5.5 Ibs. per cu.ft. V/s, 0.032 in. thick 10 pee 779 
4-B-3 5.5 Ibs. per cu.ft. 1/, 0.032 in. thick 10 ag 729 
4-C-1 5.5 lbs. per cu.ft. 1, 0.032 in. thick 15 x 583 
4-C-2 5.5 Ibs. per cu.ft. V/s, 0.032 in. thick 15 ares 548 
4-C-3 5.5 lbs. per cu.ft. 1/4 0.032 in. thick 15 er 490 
4-D-1 5.5 lbs. per cu.ft. 1/, 0.032 in. thick 20 x 554 
4-D-2 5.5 lbs. per cu.ft. VW, 0.032 in. thick 20 ee 61 
4-D-3 5.5 Ibs. per cu.ft. 1/, 0.032 in. thick 20 ae 194 
4-AA-1 5.5 Ibs. per cu.ft. 3/, 0.032 in. thick 5 x No Record 
4-AA-2 5.5 Ibs. per cu.ft. 3/4 0.032 in. thick 5 — 1,099 
4-AA-3 , 5.5 Ibs. per cu.ft. 3/5 0.032 in. thick 5 ee 660 
4-BB-i 5.5 lbs. per cu.ft.” 3/s 0.032 in. thick 10 x 610 
4-BB-2 5.5 Ibs. per cu.ft. 3/5 0.032 in. thick 10 wate 956 
4-BB-3 5.5 lbs. per cu.ft. 3/5 0.032 in. thick 10 Le 1,058 
4-CC-1 5.5 Ibs. per cu.ft. 3/5 0.032 in. thick 15 x 786 
4-CC-2 5.5 lbs. per cu.ft. 3/5 0.032 in. thick 15 ae 825 
4-CC-3 5.5 Ibs. per cu.ft. 3/5 0.032 in. thick 15 RPE Below 250 
4-DD-1 5.5 lbs. per cu.ft. 3/s 0.032 in. thick 20 x 711 
4-DD-2 5.5 lbs. per cu.ft. 3/, 0.032 in. thick 20 aaa 468 
5-A-1 5.5 lbs. per cu.ft. 1/, 52S-!/.H 5 x 1,730 
5-A-2 5.5 Ibs. per cu.ft. 1/, 0.081 in. thick 5 x 1,184 
5-A-3 5.5 lbs. per cu.ft. 1/, 0.081 in. thick 5 saa 1,435 
5-AA-1 5.5 Ibs. per cu.ft. 3/5 0.081 in. thick 5 x 1,765 
5-AA-2 5.5 Ibs. per cu.ft. 3/, 0.081 in. thick 5 x 002 
5-AA-3 5.5 Ibs. per cu.ft. 3/s 0.081 in. thick 5 63 1,775 
5-B-1 5.5 Ibs. per cu.ft, /, 0.081 in. thick 10 x 1,120 
5-B-2 5.5 Ibs. per cu.ft. /, 0.081 in. thick 10 z= 1,320 
5-B-3 5.5 Ibs. per cu.ft. 1/, 0.081 in. thick 10 re 1,420 
5-BB-1 5.5 Ibs. per cu.ft. 3/5 0.081 in. thick 10 x 1,255 
5-BB-2 5.5 Ibs. per cu.ft. 3/s 0.081 in. thick 10 x 2,055 
5-BB-3 5.5 lbs. per cu.ft. 3/5 0.081 in. thick 10 bok 1,965 
5-C-1 5.5 Ibs. per cu.ft. 1/, 0.081 in. thick 15 x 1,088 
§-C-2 5.5 Ibs. per cu.ft. 1/, 0.081 in. thick 15 ae 640 
5-C-3 5.5 Ibs. per cu.ft. 1/, 0.081 in. thick 15 aie 936 
5-CC-1 5.5 lbs. per cu.ft. 3/5 0.081 in. thick 15 x 1,416 
5-CC-2 5.5 Ibs. per cu.ft. 3/5 0.081 in. thick 15 5 1,398 
5-CC-3 5.5 Ibs. per cu.ft. 3/, 0.081 in. thick 15 bios 1,648 
5-D-1 5.5 Ibs. per cu.ft. 1/, 0.081 in. thick 20 x 975 
5-D-2 5.5 Ibs. per cu.ft. 1/, 0.081 in. thick 20 Re 872 
5-D-3 5.5 Ibs. per cu.ft. 1/, 0.081 in. thick 20 see 1,262 
5-DD-1 5.5 Ibs. per cu.ft. 3/5 0.081 in. thick 20 x 1,495 
5-DD-2 5.5 Ibs. per cu.ft. 3/5 0.081 in. thick 20 ree 1,700 
5-DD-3 5.5 Ibs. per cu.ft. 3/5 0.081 in. thick 20 : 1,182 
L 
6(U + V) = — + w" + Gb{ [2u/(b + t)] — [28u/(b + + 
0 - 
Gb{ [2u’/(b + t)] — + 2x (cof [2u/(b + t)] — + + iu) ) =0 (29) 


With 


= bw = 0 


when x = 0, one obtains after one more integration by parts 
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V) = fa ( —2Etu"iu + 2(ED), + + Gb} [2u/(b + t)] — v'} [26u/(b + + 


Gb{ [2u’/(b + — + 2x 


+ + Gb{ [2u/(b + t)] — v'} bv + + iw) ) = 0 (30) 
In Eq. (30) 6u, 6v, and 6w are arbitrary functions of the coordinate x. The equation is satisfied identically if 
—2Ejtu" + [2Gb/(b + t)]{[2u/(b + t)] — v'} = 0 (31) 
2(ED + Gb{ [2u’/(b + t)] — = 0 (32) 
2(EI) + 2x(E,/b)w = 0 (33) 
and 

2(ED + Gb{ [2u/(o + — + W = 0 (34) 

when x = L. 
2(ED ++ WwW =0 (35) 


when x = L. 

It is easily seen that the simultaneous differential equations [Eqs. (31) and (32)], together with condition (34) 
which prescribes the distribution of the external load at the free end of the cantilever, define the extensional and 
shearing deformations of the sandwich beam. The flattening of the core is governed independently by Eq. (33), 
together with condition (35). 


Solution of the Equations of Extensional and Shearing Deformations 
The general solution of the homogeneous differential Eqs. (31) and (32) is 


u =, (A;/p*) sinh px + (A2/p*) cosh px + Agx® + Age + As (36) 
b+t 1 — 1—pR\ , x8 
( rt) v= CF cosh px + A, (=) sinh px + As; (= _ 2k*x) + As (=) + Asx + Ae 
(37) 


where 
p? = (Gb)/(ED*, = 
(EI)* = + (ED, (38) 
(ED, = (1/2)t(b +E, (ED), = (1/6)fE, 


The constants of integration can be determined from the condition stated in Eq. (34) and the following end 
conditions: 


u=ov=p’' = 0, when x = 0 


L 


(39) 


= 
S 
> 
i=} 


Consideration of Eqs. (39) yields the solution 
“= (A,/p*) {sinh px — coth pL[cosh px — (p*x?/2) — 1] — Lp*x coth pL} (40) 
v = (A,/p*)[2/(b + £)]{(1 — pk*)(cosh px — 1) — (coth pL)[(1 — p*k*)(sinh px — px) — (1/6)(px)*] — 
(pL/2)(px)* coth pL} (41) 
The constant A, must be so determined as to satisfy Eq. (34). One obtains 
A; = + t)/L)pL tanh pL (42) 


where 


I, = +21; (42a) 


g 
| 
4 
— 


718 JOURNAL OF THE AERONAUTICAL SCIENCES—DECEMBER, 1948 


After substitutions and algebraic transformations, the expression for the deflection d; = v, at the free end of 
the cantilever can be written in the form given in Eq. (1). 


Solution of the Equations of the Flattening of the Core 
Eq. (33) can be written in the form 


w'Y + 4D‘w = 0 (43) 
where 
D* = (3x/bt*)(E,./E;) (48a) 
The general solution of this homogeneous differential equation is 
= B, cos Dx cosh Dx + Bs cos Dx sinh Dx + B3 sin Dx cosh Dx + By, sin Dx sinh Dx (44) 
The constants of integration can be determined from Eq. (35) and the following end conditions: 
w=w’ = 0, when x = 0 ; 
w” = 0, when x = L | 
Consideration of Eqs. (45) yields 
w = B.[cos Dx sinh Dx — sin Dx cosh Dx + (tan DL + tanh DL) sin Dx sinh Dx} (46) * 
~ (35) is then satisfied if 
= (W/ xE.)(b/L) [(DL)(cos DL) (cosh DL)/(cos* DL + cosh? DL)} (47) 


As the shortening d; of the distance between the two faces at the free end of the cantilever is twice the value of 
w at the same section, d; can be given in the form presented in Eq. (9). 


Equations of Equilibrium of a Sandwich Column 
When o + 0, g = 0, and the two ends of the column are simply supported, ep. (29a) and (29b) must be re. 
placed by the following end conditions: 
= = jw = 0 (48a) 
when x = 0, and 


= w" = iv = bw = 0 (48b) 


when x = L. Consequently, after the first integration by parts, the variation of the total potential will not 
include the last expression in parentheses in Eq. (29). On the other hand, the following term must be added to 
the right-hand member of Eq. (29) because of the variation of the potential of the end load 2ot: 


L L 
= f [(v’ + w’)(6v’ + dw’) + (v’ — w’)(6v’ — dw’)dx] = (v"év + w"dw)dx (48c) 
0 


After the second integration by parts, the right-hand member of Eq. (30) must be replaced by its first term in 
parentheses augmented by the right-hand member of Eq. (48c). Hence, the equations of equilibrium become: 


—2E,tu” + [2Gb/(b + #)]{[2u/(b + t)] — v’} = 0 (49) 
+ — Gb)v” + Gb[2u’/(b + = 0 (50) 
2(EI) + 2otw"” + 2x(E,/b)w = 0 (51) 


“Again u and v are interconnected by the simultaneous linear differential equations [Eqs. (49) and (50)], while 
w can be calculated independently from Eq. (51). 


Buckling Without Flattening 
After algebraic transformations Eqs. (49) and (50) can be written in the form 


— + [Gb/(ED)Jy’ = 0 (52) 
+ [(2ot — Gb)/2(ET) + [Gb/2(ED),\u’ = 0 (53) 
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y = (1/2)(6 + (53a) 
The general solution of these equations is 
u = A, cosh pix + Az sinh pix + As cos pox + Ag sin pox + As (54) 
y= sinh + A, cosh + + as sin pox — A4gcos pox) + Asx + Ace (55) 
where 
pr =|Va—>| 
= (1/2)| Ve? » = (1/2)8 
8 = [(Gb)(EI), — = [(Gb)(20t)]/ (ED, (2ED);) 
Because of the simple supports the end conditions are: p? = (1/2)[|WN| — (4/2)] (63a) 
tes = (1/2)[|VN| + (M/2)] 
when x = 0 and when x = L. Consequently, provided 
=y=y" =0 (57b) D<0O (63b) 


when x = 0 and when x = L. These end conditions w = 4, cos px + Az sin px + Asx cos px + Aux sin px 
can be satisfied by Eqs. (54) and (55) only when all 


the integration constants vanish (trivial solution) or (64) 
when where 
(58) p? = (M/2) =|VN| (64a) 
A i = 
provided 
The nontrivial solution exists, therefore, when 
D=0 (64b) 
pol => nT (59) 
and 
where ” is any integer. Substitutions from Eqs. (56), 
transformations, and solution for P,, = 2et yield w = A, cos px + Az sin px + A; cos gx + Ag sin gx 
(65) 
2 
n?P, + Gb where 
When the two ends of the column are pinned, m = 1, Pp? = (1/2)[M + (M? — 4N id (65a) 
and the buckling load is given by Eq. (12). The 
symbols P; and P2 are defined in Eqs. (12a) and (12b). q’ = (1/2)[M (M 4N)") 
provided 
Buckling with Flattening 
Eq. (51) can be rewritten in the form 
With the two ends of the column simply supported, 
w + Mw + Nw 0 (61) the end conditions are 
where w=w" =0 (66) 


M = ot/(EI), N = (x/bI;)(Ee/Ey) (61a) when x = 0, Z. These end conditions admit only of a 


The solution of this homogeneous differential equa- trivial solution when D < 0. In the other two cases 


tion depends upon the value of the discriminant sinusoidal deflections are possible when 
“D = M?—4N : (62) sin pL = 0 (67) 
The solution is sin gl = 0 
w = A, cosh px cos gx + A: cosh px sin gx + that is, 
A; sinh px cos gx + A, sinh px sin gx (63) pL = mn (68) 
where gL = nx 
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Fic. 20. Constants needed in calculating buckling with flatten- 
ing. 


Substitution of the values of p and gq from Eqs. (65a) 
yields the buckling condition 

M = (n/L’)? + (L'/m)?N (69) 
where 


L' = L/n (69a) 


is the wave length at buckling. Of all possible wave 
lengths the one is of practical interest that makes 
the buckling stress—which is proportional to \/—a 
minimum. The condition of this minimum is 
—2x? 2L’ 
aL’, 
After algebraic transformations, Eq. (70) can be 
written in the form 
(3b/2L')? = 


(70) 


(71) 
where 


= (3/4)*[x + (71a) 


and the dot over x signifies differentiation with respect 
to (2L’/rb). 
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Values of ® are plotted against (7b/2L’)? in Fig. 20. 
With the aid of this figure it is easy to find that value 
of (rb/2L’)? which satisfies Eq. (71). For con- 
venience, Eq. (69) can be rewritten in the following 
form: 


1/t\?/ rb \? 1/b\ /2L’\? 
From Eq. (72) the critical stress can be easily calculated 
with the aid of the value of (1b/2L’)? just determined 


and the values of x. which are also plotted against 
(rb/2L’)* in Fig. 20. 


In the case when D = 0, the solution reduces to 


Ser = (73) 


The writers wish to mention that, although it is 


presented in a different form, the solution of the prob- 


lem of buckling with flattening derived here is the same 
as that given by them in reference 2. 
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N.A.C.A. Investigations of Gas-Turbine- 
Blade Cooling 


HERMAN H. ELLERBROCK, JR.* 


National Advisory Committee for Aeronautics 


ABSTRACT 


The effect of increasing turbine-inlet temperature on turbine- 
propeller- and turbojet-engine performance is discussed and 
marked improvements are noted. Limitations to the increase in 
turbine-inlet temperature of present-day engines are brought out, 
and the need for turbine-blade cooling and the status of this 
problem are discussed. 

Analytical studies on turbine-blade cooling made by the 
N.A.C.A. are discussed, and brief descriptions of results 
obtained in the studies comparing the increase in effective gas 
temperature obtainable with different methods of cooling are 
given. Because of the difficulty, if not impossibility, of solving 
turbine-cooling problems theoretically, recourse is had to experi- 
mental investigations. A general discussion of the use of rigs of 
static cascades and rotating cascades of blades to obtain experi- 
mental data is given. 


SYMBOLS 


a = speed of sound in air, ft. per sec. 

A = area of blade cross section, sq.ft. 

Cys = specific heat of fluid, B.t.u. per (Ib.)(°F.) 

d = hydraulic diameter, ft. 

dA, = infinitely small area of blade surface, sq.ft. 

dAw = infinitely small area of turbine casing wall, sq.ft. 


g = acceleration of gravity, ft. per sec.” 

Gr = Grashof Number, gd*8py?(t — ty) /uy? 

he = convection heat-transfer coefficient, B.t.u. per (sq.ft.) 
( °F.) (sec.) 


he. 9 = convection heat-transfer coefficient from gas to blade, 
B.t.u. per (sq.ft.)(°F.) (sec.) 

ke = thermal conductivity of ceramic coating on outside 
blade surface, B.t.u. per (hour) (sq.ft.)(°F. per ft.) 

ky = thermal conductivity of fluid, B.t.u. per (sec.)(sq.ft.) 
(°F. per ft.) 

km |= thermal conductivity of blade metal, B.t.u. per (sec.) 
(sq.ft.)(°F. per ft.) or B.t.u. per (hour)(sq.ft.) 
(°F. per ft.) 

l = length of any passage, ft. 

i = blade length from root to tip, ft. 

M = Mach Number, V;/a 


Ne = polytropic exponent of compression 

Nu = Nusselt Number, h.d/ky 

p = peritneter of blade section, ft. 

P, = total pressure at compressor inlet, lbs. per sq.ft. 

P, = total pressure at compressor outlet, lbs. per sq.ft. 

P, — = total pressure at turbine outlet, lbs. per sq.ft. 

Pr = Prandtl Number, ¢p,su;/ky 

q = net heat transfer to blade from wall radiation for small 


finite area, B.t.u. per sec. 
distance between points for which net radiation is 
being determined, ft. 


ll 
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Re = Reynolds Number, Vypsd/uy 

t = temperature of any heated or cooled section, °F. 

ly = mean fluid temperature, °F. 

tn = blade-metal temperature, °F. j 

tm.r = blade-root temperature, °F. 

t.e = gas temperature affecting heat transfer (effective gas | 
temperature), °F. 

At,.¢ = increment of effective gas temperature, °F. 

T, = temperature of finite area of blade receiving radiant 
heat from wall, °R. 

Tw. = temperature of finite area of turbine-casing wall radiat- 
ing heat to blade surface, °R. 

T, = total temperature of gas at compressor inlet, °R. 

T; = total temperature of gas at turbine inlet, °R. 

T,.e« = gas temperature affecting heat transfer, °R. 

T,.2 = static temperature of gas, °R. 

T,;, = total temperature of gas, °R. 

V = rotational speed of blade, ft. per sec. 

Vmar. = Maximum rotational speed at which blade can be 
safely operated, ft. per sec. i 

V; = velocity of fluid, ft. per sec. f 

x = distance from blade root to point considered on blade, 


ft. 
Xe = critical distance from blade root to point considered on 
blade, ft. 
a = factor relating effective, static, and total gas tempera- 


tures 

a, = absorptivity of turbine-blade surface : 

a@w = absorptivity of turbine-casing wall 

B = coefficient of thermal expansion of fluid, (°F.)~! 

nm,¢ = mechanical efficiency of compressor 

nm.t = mechanical efficiency of turbine 

uy = viscosity of fluid, slugs per (ft.) (sec.) 4 

py = density of fluid, slugs per cu.ft. 

o = Stefan-Boltzmann constant for radiation, B.t.u. per 
(sec.) (sq.ft.)(°R.)4 

¢» = angle between normal to infinitely small area of blade 
surface and line connecting this area with infinitely 
small area of turbine-casing wall for which radiant 
heat transfer is evaluated, deg. 

¢%» = angle between normal to infinitely small area of 
turbine-casing wall and line connecting this area 
with infinitely small area of blade surface for which 
radiant heat transfer is evaluated, deg. 


INTRODUCTION 


—_ ESSENTIAL OBJECTIVE of turbine research in the 
; field of aeronautics is the improvement of pay 
loads, range, and reliability of aircraft. The greatest 
efforts in turbine research have therefore been devoted 
to increasing the power that can be developed by 
an engine of a given size and weight, to increasing 
the power that can be realized from a given quan- 
tity of fuel, and to improving the reliability of tur- 
bines. 


20. 
ue 
ng | 
72) 
ted 
ied 
nst 
73) 
is 
ob- 
me 
‘aw: 
3, 
[olt, 
num 
any 
ung, 
nter- 
‘eins 
S., 
any, 
bical 
num 
. 15, 
|_| 


Although substantial increases in the power and the 
economy of turbines can be realized from improvements 
in the flow capacity and the efficiency, the largest and 
most significant benefits may be derived from increas- 
ing the temperature of the working fluid at the turbine 
inlet. For a given capacity or airflow quantity, it can 
be shown from turbine theory that the power per pound 
of air is proportional to the turbine-inlet temperature. 

Present hydrocarbon fuels can produce far higher 
temperatures than are now used, but current turbine 
materials have insufficient strength at these higher 
temperatures. While better heat-resisting materials 
are being developed, turbine-blade cooling is a more 
immédiate means of realizing the remarkable perfor- 
mance improvement that should be possible with 
raised turbine-inlet temperatures. 

Because of this indicated large improvement in per- 
formance due to increased turbine-inlet temperature, 
the National Advisory Committee for Aeronautics 
started investigations on cooling of turbine blades at its 
Cleveland laboratory in 1945 to obtain ultimately tur- 
bines that will operate at temperatures much higher 
than those allowable at present. 

This paper briefly describes some aspects of the 
N.A.C.A. investigations. The effect and importance 
of increased turbine-inlet temperature and the need for 
turbine cooling are demonstrated. The status of tur- 
bine-cooling research and some of the problems inherent 
in the investigations are mentioned to indicate the 
amount of research effort necessary. Analytical inves- 
tigations and their results are given, and the experi- 
mental procedures needed to supplement them are 
described. 

The investigations reported are, for the main part, 
the result of work by the following members and former 
members of the Cleveland laboratory staff: Dr. W. 
Byron Brown, Dr. John N. B. Livingood, Mr. William 
D. Bowman, and Mr. Lincoln Wolfenstein. 


EFFECT OF INCREASED TURBINE-INLET TEMPERATURE 
AND NEED FOR TURBINE-BLADE COOLING 


Thermodynamic-cycle analysis of the simple turbojet 
propulsion system indicates that, if the turbine-inlet 
temperature is raised to 2,000°F. from a temperature of 
1,500°F., the specific thrust of the engine may increase 
32 per cent, although at the expense of a 7 per cent 
increase in specific fuel consumption because of loss of 
propulsive efficiency at the high relative jet velocities. 


At the same time, the specific weight of the power’ 


plant would be reduced approximately 30 per cent. 

For the turbine-propeller type of engine at optimum- 
power compressor pressure ratio, an increase in turbine- 
inlet temperature from 2,000° to 3,000°F. would result 
in a 57 per cent increase in specific power output with a 
5 per cent reduction in specific fuel consumption. 
Furthermore, comparison of the 3,000°F. turbine- 
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losses. Airplane speed, 
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TURBINE-INLET TEMPERATURE, °F 


Fic. 1. Turbine-propeller-engine performance without cooling 
500 m.p.h.; Mach Number, 0.69; alti- 


tude, 30,000 ft. 
propeller power plant having an optimum-power com- 
pressor pressure ratio with present-day designs shows a 
reduction of 50 per cent in specific fuel consumption and 
an increase of over 200 per cent in the specific output. 
Full utilization of high operational turbine-inlet 
temperatures, however, demands an increase in com- 
pressor and turbine pressure ratios. In the turbine- 
propeller engine, the theoretical compressor pressure 
ratio required for optimum power for a given ratio of 
turbine-inlet temperature to compressor-inlet tempera- 
ture is given by the equation: 


Py ne/{2(ne — 1)] 


A typical set of curves illustrating the effect of tur- 
bine-inlet temperature on uncooled-turbine-propeller- 
engine performance is given in Fig. 1. The specific 
thrust power and the specific fuel consumption for a 
speed of 500 m.p.h. and at an altitude of 30,000 ft. have 
been plotted for a range of turbine-inlet temperature 
from 2,000° to 4,000°F. The optimum compressor 
pressure ratio has been used for each temperature ratio. 
The rapid increase in specific thrust output with a 
slight decrease in specific fuel consumption is evident. 
Comparison of the approximate power of present-day 
engines, as indicated by the point on the figure at 1,500° 
F., with the power that could be obtained at 3,000°F. 
indicated on the curve results in a ratio of the powers of 
more than 3 to 1. Similar calculations for a turbojet 
show the same trend for specific power, but specific fuel 
consumption increases above 2,000°F. 

Although the specific power is increased appreciably 
as the turbine-inlet temperature is increased, if a turbine 
is of given design, the flow in pounds of air per second 
decreases as the turbine-inlet temperature is increased. 
The theoretical net result of changing turbine-inlet 
temperature on a turbine of a given design, if no other 
factors are changed, is that the power delivered by the 
turbine will be unchanged. Both turbine speed and 
compressor pressure ratio must therefore be increased in 
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order to retain constant flow capacity and flow simi- 
larity; this results in increasing turbine rotor and 
blade stress in order to realize the potential improve- 
ments in specific power. The increase in turbine speed 
can be minimized by altering ‘turbine geometry and 
blade configuration. 

There is every reason, on the basis of the foregoing 
discussion, for developing gas turbines capable of operat- 
ing at elevated temperatures. Present hydrocarbon 
fuels can readily produce gas temperatures in the region 
of 3,500°F. The chief obstacle to use of high tempera- 
tures, however, is that tried and tested materials having 
more strength at these elevated temperatures than 
present materials have at ordinary temperatures are 
required and are not available. Modern high-tempera- 
ture alloys lose strength rapidly at temperatures in the 
vicinity of 1,600° to 1,700°F. Moreover, the melting 
point decreases with increase in alloy content. Thus, 
the melting-point range (2,300°-2,400°F.) approaches 
the desired operating temperature as alloy content is 
increased, and the opportunity for further increase in 
strength at elevated temperatures is limited. The use- 
fulness of ceramics and other ultrahigh-temperature 
materials is still to be determined. 

Because the development of heat-resistant materials 
will require time and successful development cannot be 
predicted with certainty, the use of turbine-blade 
cooling represents a more immediate and positive solu- 
tion of the high-temperature-turbine problem and the 
realization of the remarkable performance improvement 
that has been indicated. Another important possi- 
bility of the use of blade cooling which should be men- 
tioned is that nonstrategic metals can be used for the 
turbine blades and wheels. The turbine-cooling idea is 
fundamentally the same as the cooling of those struc- 
tural parts of the reciprocating engine that are exposed 
to hot gases. In this engine, as much as 40 per cent of 
the thermal energy of the fuel has had to be dissipated 
to the cooling fluid in order to achieve the potential 
thermal efficiencies and power of the cycle. 


STATUS OF TURBINE-BLADE-COOLING PROBLEM 


At the beginning (1945) of N.A.C.A. research in 
turbine-blade cooling, no theoretical or experimental 
work was available to use as a guide. The factors that 
required study had to be determined even before the 
efficacy of various cooling methods could be worked out. 


Basic Problems 


The factors that were found to require study were 
those needed to determine some of the performance 
aspects of a cooled-turbine power plant, such as the 
losses caused by cooling, the net power output, the 
quantity of coolant required for adequate cooling, etc. 
Design conditions can also be obtained from these fac- 
tors—for example, the speed at which the engine can 
be operated satisfactorily, the gas temperature that can 


be used, the limiting length of the blade trailing-edge 
section, and the limiting length of the blade critical 
section. 

The required factors can be enumerated as follows: 


(a) Convection heat-transfer coefficient from blade 
to coolant. 

(b) Convection heat-transfer coefficient from gas 
to blade. 

(c) Radiation heat-transfer coefficient for heat 
radiated to outside blade surface. 

(d) Radiation heat-transfer coefficient from blade 
to insert inside blade. 

(e) Gas temperature affecting heat transfer (effec- 
tive gas temperature). 

(f) Blade-metal temperature. 

(g) Coolant temperature affecting heat transfer 
(effective coolant temperature). 


All problems of convection heat transfer are problems 
of flow of heat through the boundary layer. Conse- 
quently, it would be supposed that setting up the equa- 
tions for the flow of fluid and the temperature variation 
in the boundary layer could establish the heat-transfer- 
coefficient laws. Dryden! has set up the equations for 
both laminar and turbulent flows. For laminar flow, 
nine equations, five of which are differential equations, 
result with nine unknown terms. Up to the present, 
little progress has been made in the solution of this 
formidable problem. Simple problems, such as for 
flow through a pipe and across flat plates, with the 
assumption that the flow is one-dimensional or two- 
dimensional and many other assumptions, have given 
rise to a solution of the equations that fits experimental 
results fairly closely. For turbulent flow, the equa- 
tions are more difficult, the fluctuating-velocity terms 
and their equations entering into the problem in addi- 
tion to equations of the type for laminar flow. Conse- 
quently, the theory of heat transfer in eddying flow is 
unsatisfactory at present. 

As in many other fields of technical physics where the 
differential equations governing the phenomena are 
such that their solution is impractical, if not altogether 
impossible, experimental techniques aided by the 
method of dimensional analysis based on the 7 theorem 
have proved a valuable tool in the treatment of heat- 
transfer problems. By use of such an analysis, the 
Nusselt Number, which involves the convection heat- 
transfer coefficient, can be shown to be a function of 
the parameters given in the following equation: 


Nu=y Re, Pr, Gr, M, B(t — “| (2) 
My Uy 


For general forced-convection heat-transfer problems 
encountered in the past, the foregoing equation reduces 
to 


"Nu = f(Re, Pr) (3) 


the other terms having little or no effect. 
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For cooled turbines, many factors may invalidate to 
some extent the heat-transfer-coefficient laws already 
established. One factor is that the flow through the 
turbine is thought to be different from the flow for which 
the laws were established. Obviously, the flow condi- 
tions in a turbine are different from those in a static rig. 


In any rig the wakes from the nozzles influence condi- * 


tions downstream; in the rotating rig an’ additional 
effect due to relative motion between nozzles and tur- 
bine blades is encountered. Radial pressure gradients 
result from turbine rotation as well as from the turning 
action of the blades. Because of the high rotational 
speeds, centrifugal forces on the gas particles are thou- 
sands of times greater than the normal force of gravity. 
Radial pressure gradients are therefore more pro- 
nounced, and cross flow is greater than in a static rig. 
Centrifugal forces affect the cooling-passage boundary 
layer as well as the gases on the outer surface of the 
blade. These phenomena have been discussed by 
Weske.? Another effect that has been generally neg- 
lected but is of particular importance in the turbine is 
the Grashof Number effect, which evaluates natural 
convection. For the turbine, the gravity term for the 
coolant can become about 10,000 to 50,000 times that 
of the earth’s force of gravity because of the rotation of 
the turbine. Consequently, the Grashof Number 
value for natural-convection circulation assumes propor- 
tions never before experienced and will probably be 
included in any equation for the inside convection heat- 
transfer coefficient of turbine blades. Radiation also 
plays a prominent part in the transfer of heat in a tur- 
bine. A difficult effect to determine and for which no 
data exist is the geometry factor that enters the 
radiation-coefficient equation. 


Evidence exists that, as the gas proceeds along the 
blade, the temperature of the gas adjacent to the blade 
surface, which affects the heat transfer, decreases be- 
cause of the action of the coolant on the blade. Thus 
the laws of the variation of gas temperature around the 
blade must be established. Recourse may be had to 
boundary-layer theory for this problem, but again the 
solutions are difficult and have to be established. 
Experimental investigations are about the quickest 
solution. 


This paper cannot give all the problems arising in the 
study of turbine cooling that are preventing the rigorous 
design of cooled turbines. Only a few have been men- 
tioned to indicate that the problem of blade cooling re- 
quires much research effort. 


Agreement of German Work 


At the end of World War II various references became 
available to the N.A.C.A. on work the Germans had 
been conducting at the same time the N.A.C.A. was 
pursuing its research. An excellent summary of the 
German theoretical investigations, which in general 
show the comparison of various cooling methods, was: 


compiled by K. J. Miiller for the Deutsche Luftfahrt- 
forschung. The agreement between these results and 
the N.A.C.A. analytical investigations is surprising, 
inasmuch as the two groups worked entirely independ- 
ently. 

It was learned from the references that the practice 
of blowing air through hollow blades was in general use 
on the latest types of German engine. Temperatures 
of about 1,550° to 1,650°F. could be attained there- 
with, even though the heat-resisting quality of the 
material used was not equal to that of American turbine 
materials. Also considered, and tested in a static cas- 
cade of blades, was the use of a thin layer of cold air 
flowing out of slits in the blade and along the outside 
blade surface. This method is called ‘‘film cooling.” 
This protective layer acted as a shield between the hot 


_ gases and the blade surface. A favorable cooling effect 


with a small amount of air was obtained. The greatest 
water-cooling progress was made with a single-stage 
turbine that ran in continuous operation at a turbine- 
inlet temperature of 2,200°F. with a cooling process 
devised by Dr. Ernst Schmidt of the Hermann Goering 
Institute. 


ANALYTICAL INVESTIGATIONS 


Many methods have been suggested for accomplish- 
ing the desired cooling of turbine blades. The decision 
as to which of these methods or what combination of 
these methods has the best practical value can be made 
only after a large amount of experimental data has been 
obtained and studied. The mechanism of flow in 
boundary layers surrounding turbine blades is not yet 
sufficiently understood to make accurate quantitative 
comparisons by analysis. If the cooling of the bound- 
ary layer is neglected, however, a great deal can be 
learned about the suitability of various forms of conduc- 
tion cooling by employing the principles originally set 
forth by Fourier. Suitable average values may be 
assigned to the heat-transfer coefficients between the hot 
gases and the metal and between the metal and the 
cooling fluid. Allowances may also be made for the 
effects of radiation. 

The N.A.C.A. applied this procedure to three possible 
methods of conduction cooling: (1) rim cooling, in 
which the temperature of the rim of the turbine disc is 
reduced by the cooling medium and the solid blades are 
cooled by conduction of heat to the rim; (2) hollow- 
blade air cooling, in which air is passed through hollow 
turbine blades and heat flows through the thin metal 
shell forming the turbine-blade contours to the air; 
and (3) liquid cooling, in which liquid is circulated 
through passages formed in the turbine blades and discs 
and heat flows through the metal to the cooling fluid. 


Method of Analysis 


In order to evaluate temperature distribution in a 
cooled turbine blade, the following parameters must be 
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specified: the turbine structures; the shape and the 
material of the blades; and the temperatures, the 
rates of flow, and the composition of the hot gases and 
the cooling fluid. The heat-transfer and radiation 
coefficients can then be determined if their laws are 
known. Because of lack of data, the equation for fully 
developed turbulent flow inside pipes has been used to 
determine an average convection heat-transfer coeffici- 
ent from metal to coolant, and, in general, the equation 
for flow across single cylinders @vith the outside blade 
perimeter divided by 7 as the dimension in the Reynolds 
and Nusselt Numbers) has been used to determine an 
average convection heat-transfer coefficient from gas to 
metal. The usual equations for radiation were used to 
evaluate these effects with a geometry factor of 1.0. 
Aljl these formulas were obtained from McAdams’ heat- 
transfer text. The temperature of the gas affecting 
heat transfer was assumed to be constant around the 
blade and was determined from the equation: 


(Ty,¢ Ty,s)/(T =) a8) (4) 


where a is roughly equal to the square root of the 
Prandtl Number of the gas.‘ The static temperature 
is determined from the turbine-inlet total temperature, 
the gas-flow rate, and the turbine-nozzle area. The 
total temperature 7,, is determined from the static 
temperature and the relative velocity of the gas to the 
blade. - 


From the foregoing data, the temperature distribu- 
tion for any cooling arrangement and at any gas tem- 
perature can be found by formulating the equations 
expressing the heat balance at each infinitesimal part of 
the turbine and performing the required mathematics. 
In order to determine the amount by which the tempera- 
ture of thc hot gases can be increased by a given method 
of cooling, the temperature distribution along the blade 
radius must be found and compared with the allowable 
temperature distribution determined from the stress in 
the blades and stress-rupture data. The desired life of 
the blade is then assumed, and the relation between the 
maximum allowable temperature of the blade material 
and the blade stress is established as follows: When the 
tip speed of the blade is specified, the radial stress dis- 
tribution of the blade is found, and the maximum per- 
missible temperature at each radius of the blade may 
then be computed. A curve may then be drawn to 
show the radial distribution of the stress-limited tem- 
perature. 


Representative curves of the stress-limited tempera- 
ture distribution for several tip speeds are shown in 
Fig. 2, together with a typical curve of the computed 
temperature distribution. The solid curve shows the 
temperature along the blade as determined by heat- 
transfer equations. The dotted curves show the 
allowable stress-limited temperatures. For simplicity, 
the stress distribution used in this example is that due 
only to centrifugal stress in a uniform tapered blade. 
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Fic. 2. Method of determining limiting speed and critical 
blade section from curves of temperature distribution and allow- 
able stress-limited blade temperature. 


The various other important stresses, such as bending, 
vibratory, and thermal, have been neglected; thus the 
life of the blade is assumed to be limited only by the 
stress-rupture characteristics of the material. The 
maximum gas temperature for a given tip speed is estab- 
lished when the temperature-distribution curve has a 
point of tangency with the stress-limited distribution 
curves. The point of tangency is called the critical 
point. When the curves intersect, the sections on the 
blade within the overlapping regions are too highly 
stressed. The effectiveness of a given method of cool- 
ing is determined by the increase in gas temperature 
produced by cooling for a given stress-limited distribu- 
tion curve. The curves may also be used to determine 
the maximum tip speed for a given calculated tempera- 
ture distribution. 


Rim Cooling 


The effectiveness of rim cooling—that is, cooling the 
rim of the disc so that the heat flows through the solid 
blades to the rim—was investigated for a wide range of 
cooling variables and blade shapes, and the results are 
summarized in Fig. 3. The allowable increase in the 
effective gas temperature is plotted as ordinate against 
the difference between the effective gas temperature 
and the temperature at the blade root for each of a series 
of cooling patameters. The dimensionless parameter 
includes gas-to-metal heat-transfer coefficient /,,,, the 
perimeter of the blade section p, the blade length L, 
the thermal conductivity of the blade material k,,, and 
the cross-sectional area of the blade section A. For 
modern gas turbines of high power and heat-transfer 
coefficient, the perimeter and the blade length are 
large, whereas the thermal conductivity is small. For 
example, the value of the parameter would be 14 for 
turbines similar to one of our modern jet engines. For 
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able Mach Number, 0.5. 


turbines similar to the B-2 turbosupercharger, the value 
is about 4.5. 

Under these circumstances, the cooling is indicated 
by the lower curves, and the allowable increase in gas 
temperature by rim cooling rarely exceeds 200°F. 
Although this temperature increasé is sufficient to cause 
a significant gain in turbine performance, the increase 
is not nearly enough to permit use of the fuel-limited 
temperature of about 3,500°F. By the use of short, 
thick blades of high-conductivity materia) (which is not 
available at present with the desired stress-rupture 
characteristics), the effectiveness of rim cooling would 
be indicated by the upper series of curves. Even under 
these extremely favorable circumstances an effective- 
gas-temperature increase of only about 675°F. could be 
obtained. The study also showed, however, that, with 
uniform strength characteristics, blade life could be 
increased appreciably by small decreases in rim tem- 
perature at a fixed gas temperature. 


Air Cooling 


The effectiveness of the more direct method of using 
air by passing it through the center of a hollow blade 
similar in design to that shown in Fig. 4 is shown by 
the lower curve in that figure. The permissible increase 
in effective gas temperature is plotted as ordinate 
against dilution (the ratio of the cooling-air flow to the 
hot-gas flow). The curve terminates at the point 
where choking occurs in the coolant passage or the air 
reaches a Mach Number of 1.0. When the cooling-air 
Mach Number was limited to 0.6, an increase in effec- 
tive gas temperature of about 500°F. was found for the 
hollow blade. The cooling-air quantity required for 
this increase was abdut 13 per cent of the hot-gas flow. 
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The analyses showed that the critical section of the 
blade occurs close to the root section. Only the cooling 
air immediately adjacent to the internal blade surface 
is effectively utilized for cooling the blade in the critica] 
regions. The necessary high cooling-air velocities 
therefore might possibly be maintained, and the airflow 
quantity required might be reduced by using an insert 
of the type shown in Fig. 4 in the hollow blade, with air 
passing through the space between the insert and the 
blade. The results of a study with this method of 
cooling are also shown in Fig. 4. The flow of cooling 
air for the same critical blade conditions was found to be 
reduced to approximately one-third of the flow without 
the insert. An increase in effective gas temperature of 
650°F. was found to be possible with a dilution of , 
slightly more than 4 per cent and a cooling-air Mach 
Number of 0.6. These increases in effective gas tem- 


_ perature are much larger than the increases of the order 


of 100°F. that would be available with rim cooling for 
similar blades or even than the 250°F. possible by the 
use of high-conductivity turbine-blade metals if they 
were available. With a cooling-air Mach Number of 
1.0 and for the blades used, the methods of air cooling 
studied to date indicate that gas temperatures up to 
2,300°F. at the turbine inlet are permissible. 


Liquid Cooling 


When liquids are used for the cooling of turbine dises 
and blades, a closed cooling system is probably neces- 
sary to avoid wasting the liquid; otherwise, liquid- 
coolant.consumption would be of the same vital interest 
as fuel consumption. For this system of cooling, the 
liquid could flow radially outward through one or more 
passages in the blade and radially inward through one 
or more. Because of the large centrifugal force, the 
pressures of the cooling liquid would be extremely high. 
In designing the cooling passages, the stresses induced 
in the blades by the coolant must therefore be con- 
sidered. A study has shown that the passages should 
have a circular cross-section to minimize bending stresses 
and that the metal surrounding the cooling sections 
should have a considerable thickness. 
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Fic. 4. Variation of permissible gas-temperature increase with 
dilution. Life, 1,000 hours; blade Mach Number, 0.5. 
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THERMAL CONDUCTIVITY, kp, , 5 BTU/(HR)(SQ FT) (°F/FT) 


THERMAL CONDUCTIVITY, km , (00 BIU/(HR) (SQ FT) (°F/ FT) 


L100 
tm. °F 


650 


Fic. 5. Isotherms in blade sections of materials with different 
conductivities with liquid cooling. Gas flow, 55 Ibs. per sec.; 
water flow, 6.42 Ibs. per sec.; gas temperature, 2,000°F.; water 
temperature, 200°F. 3 


Among the studies made on liquid cooling was that of 
determining the temperature gradients in a plane per- 
pendicular to a radius. The gradients were calculated 
by relaxation methods. Some typical results are shown 
in Fig. 5. The distributions shown are for thermal 
conductivities of 15 (representative of modern high- 
temperature alloys) and 100 (approximately the con- 
ductivity of aluminum) B.t.u. per (hour) (sq.ft.) 
(°F. per ft.). If the blade is completely surrounded by 
gases of uniform temperature, excessively high tem- 
peratures will obviously be encountered at the trailing 
edge of a blade of the configuration shown, particularly 
when the metal thermal conductivity is low. As re- 
ported previously, there is evidence that the gas tem- 
perature cools in the boundary layer as the gas pro- 
ceeds along the blade; consequently, the trailing-edge 
temperatures shown probably correspond more nearly 
to turbine-inlet temperatures of 3,000°F. rather than 
the 2,000°F. used in the analysis. Experimental data 
are needed before more accurate estimates can be made 
of boundary-layer effect. 

Further analytical studies showed that the addition 
of cooling passages in the vicinity of the critical regions 
reduced the trailing-edge temperature. On the basis of 
these studies, if the cooling passages are placed suffi- 
ciently close to the trailing edge of the blade, about '/, 
in., sufficient cooling could be obtained to permit the 
operation of modern turbine materials at least to 
turbine-inlet temperatures of 3,000°F. Analyti- 
cal investigations of the type made to obtain the results 
of Fig. 4 have to be made before this value can be veri- 
fied. 

Because of the rather high trailing-edge temperatures 
obtained with liquid cooling with some configurations 
when holes close to the trailing edge are impractical, a 
short study was made to determine the effectiveness of 


ceramic coatings on water-cooled blades. Some of the 
results are shown in Fig. 6 for the blade configuration 
illustrated. The increase in effective gas temperature 
over that obtainable with liquid cooling alone is plotted 
as ordinate against ceramic-coating thickness as abscissa. 
Calculations were made for a point near the middle of 
the trailing part of the blade (presumably nearly the 
hottest part), as indicated by the cross mark on the 
blade configuration, for both high- and low-conductivity 
ceramics and blade metals. The results showed that 
ceramic coatings of about 0.010-in. thickness and witha 
low conductivity (2.5, about equal to that of thorium 
oxide) placed around blades made of present-day tur- 
bine materials (conductivity of 15) would permit an 
increase of about 450°F. in effective gas temperature 
above that permitted by the liquid cooling alone. 
Ceramic coatings also offer the advantage of protecting 
the metal from the erosive effects of the gases. 
Although a high degree of cooling can be anticipated 
from the use of liquid cooling, objections to this method 
of cooling might be raised on the grounds that a dis- 
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Fic. 6. Variation of increase in effective gas temperature with 
ceramic-coating thickness for two metal and two ceramic thermal 
conductivities. 
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Fic. 7. Turbine-propeller-engine performance with and with- 


out cooling losses. Airplane speed, 500 m.p.h.; Mach Number, 
0.69; altitude, 30,000 ft. 
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proportionate amount of heat must be removed from 
the cycle to achieve the required degree of cooling. In 
order to determine the losses that might be expected to 
result from this abstraction of heat, an analysis has been 
made to determine how the specific output and overall] 
brake efficiency, among the many factors studied, of a 


turbine-propeller engine vary with the turbine-inlet - 


temperature when the turbine is liquid-cooled. The 
amount of heat that had to be abstracted by the cooling 
fluid to maintain reasonable temperatures in the tur- 
bine blades was calculated, and due allowances were 
made for the losses occasioned by pumping the fluid 
through the passages and the losses chargeable to the 
internal drag of the necessary coolers. Some essential 
results of this study are shown in Fig. 7, where the 
specific thrust power and specific fuel consumption are 
plotted for the same conditions as used in Fig. 1. The 
results shown in Fig. 1 have been plotted in Fig. 7 for 
comparison. In spite of the cooling losses, the power 
available from a unit mass flow of air increases steadily 
with an increase in turbine-inlet gas temperature, and 
the specific fuel consumption decreases slightly. The 
relatively low fuel consumptions indicated in the figure 
result primarily from the use of high cycle pressure 
ratios. Comparison of the results of Figs. 1 and 7 (the 
difference in the results is caused by the cooling losses) 


shows that the losses are a small percentage of the power" 


output. This situation is comparable with that of 
reciprocating engines where the loss of heat to the cool- 
ing fluid is of secondary importance in comparison with 
benefits of increased power and reliability obtainable 
by reducing the metal temperatures. 


Needed Research 


The foregoing discussion gives some of the highlights 
of the analytical studies made to date. Further studies 
of the boundary-layer phenomena, film cooling, the 
effect of placing fins in hollow blades to augment the 
cooling surface when air cooling is used, and the effects 
of cooling of various kinds on turbine-propeller- and 
turbojet-engine performance still have to be made. 
Also, the effects of blade cooling on airplane perform- 
ance must be evaluated. The course of experimental 
research can be effectively guided by means of such 
studies. 


EXPERIMENTAL PROCEDURES REQUIRED 


In the section ‘‘Basic Problems,’’ some of the un- 
known factors in the turbine-blade-cooling problem 
were discussed. In the preceding section “Analytical 
Investigations,” the analyses had to be based on several 
assumptions. Experimental procedures are necessary 
to investigate and verify these factors and assumptions, 
for the turbine-cooling problem is extremely difficult, if 
not impossible, to solve wholly theoretically. As yet, 
however, no N.A.C.A. experimental results have been 
published. 
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Static Cascades 


In experimental work, the stationary cascade is an 
essential part of a research program leading to the 
development of reliable methods and physical rela- 
tions for rational design of cooled turbine engines. An 
adequate understanding of the advantages and limita- 
tions of blade cooling requires research that cannot very 
well be conducted on the full-scale turbine. Blade- 
cooling investigations in the actual turbine suffer four 
severe limitations: 

(1) Inaccessibility for necessary instrumentation. 

(2) Impossibility of isolating and individually con- 
trolling convection, radiation, and conduction without 
exceeding construction and operating restrictions. 

(3) Restriction of maximum gas temperature in 
existing machines, which effectively prevents rapid 
exploration of cooling processes at the high temperature 
levels. 

(4) Impossibility of rapidly investigating various 
cooling methods, blade configurations, etc. 

Both basic and empirical data can be obtained with 
static cascade rigs. The distinction between basic and 
empirical data depends on the method of approach in 
layout and operation of the rig. Basic data for the 
heat-transfer laws may be obtained by breaking down 
the problem to eliminate or control each effect and thus 
determine true convection and radiation coefficients, for 
example. The analytical methods are then used to 
combine the basic elements for design purposes. The 
empirical data would be obtained where the rig was 
designed and operated to reproduce the conditions of an 
actual turbine stage. The data would consist of tem- 
perature-response relations resulting from changes in 
inlet temperature, radiation-source temperature, angle 
of attack, and so forth. 

With static cascade rigs the effects of all or some of 
the following factors on heat-transfer coefficients and 
blade temperatures can be studied: 


(1) 


Gas-temperature variation. 


(2) Coolant-temperature variation. 

(3) Coolant flow. 

(4) Gas flow. 

(5) Rim, air, liquid, film, and boundary-layer cool- 
ing. 

(6) Blade configuration. 

(7) Blade size (blades range in size from */,- to 6- 


in. span, for instance). 
(8) Blade angle of attack. 


(9) Radiation heat transfer. 
(10) Convection heat transfer. 
(11) Ceramic coatings. 
(12) Coolant Mach, Reynolds, and Prandtl Num- 
bers. 
(13) Gas Mach, Reynolds, and Prandtl Numbers. 


(14) Blade material. 


The rigs should be so constructed that either local 
heat-transfer coefficierits around the blades, average 
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coefficients, or both can be obtained. In general, much 
of the basic data should be obtained using hot and cold 
air or cold water (which should be checked at high 
temperatures with combustion gases) in order to isolate 
convection effects, to overcome the lack of knowledge 
about combustion-gas properties, and to eliminate as 
much as possible the effects of carbon deposits on the 
blades and the possibility of incompletely burned gases. 


On the coolant side, the blades should be instrumented 


in such a manner that entrance effects on the heat- 
transfer coefficient can be obtained. The coolant flow, 
of course, is not fully developed into turbulent or 
laminar flow until it proceeds some distance through the 
blade coolant passages. Boundary-layer temperatures 
and velocities are required for a full solution of the 
problem. The use of a hot-wire anemometer is sug- 
gested. 


On the basis of aerodynamic data on cascades ob- 
tained by others, an odd number of blades is considered 
better than an even number. If one blade is instru- 
mented, at least three blades should be used; if two 
blades are instrumented, five blades; etc. Heat dams 
should be provided to eliminate heat losses from the 
blade ends as much as possible. Thermocouples should 
be so located on the blades that the spanwise distribu- 
tion around the periphery is obtained, as well as a de- 
tailed peripheral distribution along the span. In addi- 
tion, the coolant passage should be fully instrumented 
along the camber line of the blade at about two stations 
representing inlet and outlet of the instrumented part 
of the blade. An additional probe should be used in 
conjunction with these surveys to determine the condi- 
tions in the passage between the blades. 


The practical size of static rigs precludes use of pres- 
sure tubes in the boundary layer; multiplication of the 
estimated boundary-layer thickness of the gas film 
around the blades by a factor’of at least 4 would permit 
such use. Near leading edges, hot-wire probes are also 
severely limited in application. Because of aerody- 
namic loading at velocities in excess of 100 ft. per sec., 
wire length-diameter ratios must be restricted to values 
of the order of 60. End-conduction effects are large 
and difficult to calculate accurately because of uncer- 
tainties in flow and temperature distributions. Other 
efiects too numerous to mention affect the accuracy of 
measurement. Some simple calculations show that 
interferometric methods can be used to obtain 
boundary-layer measurements in static rigs of the size 
that is practical. Excellent photographs of flow around 
turbine-blade profiles were obtained by Vietinghoff- 
Scheel’ with an interferometer and further indicate 
that such methods can be used to determine boundary- 
layer characteristics. 


In order to evaluate better the effects of radiation for 
a complicated structure such as a turbine wheel or cas- 
cade of blades, a mechanical integrator’ to determine 
geometry factors can be used. It can be shown that 
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The quantity under the integral sign is equal to the geo- 
metry factor times the small finite area on the blade sur- 
face being evaluated. The geometry factor is deter- 
mined graphically by use of the integrator.? In order 
to use the integrator, large models of the rig (about 20 
times size) must be made so that the integrator can be 
placed between the blades. The geometry factor and ¢ 
are found for small finite areas over the entire blade sur- 
face. All the values of g are then summed to get the 
total heat radiated. ‘ 

Before test data are obtained on static cascade rigs, 
uniform temperature, pressure, and velocity profiles 
must be ensured for both the gas and the coolant up- 
stream of the blades. Velocity, pressure, and tempera- 
ture profiles; heat losses off the ends; and other items 
must be checked to obtain good aerodynamic and heat- 
transfer data. 


Rotating Cascades 


The advantage of being able to apply data from static 
cascade rigs, which are simple to construct and operate, 
to turbine design is evident. The applicability of the 
data is determined by checking against data from rotat- 
ing cascades of blades. Rotating cascades of blades 
are expected to shed considerable light on such factors 
as boundary-layer flow with centrifugal forces and 
under radial pressure variations on both inside and out- 
side surfaces of the blades, natural convection effects. 
tending to augment the internal heat-transfer coeffi- 
cient, compressibility temperature-rise effects in cooling- 
air passages, actual losses incident to pumping air 
through blades, and relative losses of various discharge 
methods (i.e., from tip, with and without controlled 
vector, or into boundary layer). None of these vital 
factors can be handled in static cascades. 

A turbine rig, including the blades, should be thor- 
oughly instrumented with both thermocouples and 
pressure tubes. Measurements should be taken in the 
blade-cooling passages to determine flow and heat- 
transfer coefficients. The thermocouple and pressure 
leads from rotating parts will probably have to run 
radially inward between the discs into a hollow instru- 
ment shaft, which will also serve to drive pressure and 
temperature pickups. 


SUMMARY 


The work done so far by the National Advisory Com- 
mittee for Aeronautics on turbine-blade cooling has 
been both analytical and experimental. Analytical 
studies involving several assumptions show that re- 
markable improvements in gas-turbine-engine perform- 
ance can be obtained by increasing the gas tempera- 
ture at the turbine inlet even when cooling losses are 
considered. Rim cooling with present-day metals 
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allows little increase in effective gas temperature over 
that obtainable with no cooling, but blade life can be 
increased appreciably with small decreases in blade 
temperature caused by rim cooling. At a Mach 
Number of 0.6 an effective-gas-temperature increase of 
about 500°F. can be obtained with hollow-blade air 
cooling and an increase of about,650° with air cooling 
when an insert is placed in the hollow blade, liquid- 
cooled high-melting-point alloy-steel turbines appar- 
ently should be able to operate at turbine-inlet tem- 
peratures near those obtainable with present-day fuels. 

Because of the difficulty of solving the turbine-cool- 
ing problem entirely theoretically, recourse is had to 
experiments. By means of rigs of static and rotating 
cascades of blades and experimental procedures briefly 
described, data are believed to be obtainable which, in 
conjunction with the theory of turbine-blade cooling, 


will lead to methods for the rigorous design of cooled — 


turbines. 
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Fracturing Characteristics of 
Aluminum-Alloy Plate’ 


L. J. KLINGLER? G. SACHSt+ 
Case Institute of Technology , 


ABSTRACT 


Commercial hot-rolled 1!/,-in. 24ST aluminum plate exhibited 
large variations in fracture stress and reduction in area for 
different orientations of the test specimens. These variations 
were attributed to mechanical anisotropy in the plate. Measure- 
ment of the angles at which the specimens fractured indicated the 
presence of a plane of weakness, but no fracture function could 
be found to correlate this with the fracture stresses and ductilities. 


INTRODUCTION 


WO TYPES OF ANISOTROPY may be encountered in 

annealed metal. ‘Crystallographic anisotropy” 
results from a preferred orientation or “‘texture’’ in the 
crystal structure of the metal, while ‘‘mechanical 
anisotropy” results from a geometrically directional 
arrangement of certain phases, inclusions, porosities, 
cracks, etc. In addition to these two, any cold-worked 
metal develops a third type of anisotropy, which may 
be termed ‘‘anelastic anisotropy.” This type appears 
to be caused by a certain kind of residual stress set up 
during the cold working. It can be reduced and 
seemingly eliminated by a stress relief anneal. A paper 
has been published recently’ which discusses this type 
of anisotropy. 

This paper will deal with the effect of the mechanical 
anisotropy on the fracturing characteristics of hot- 
rolled aluminum plate. A previous paper? dealt with 
the effect of the crystallographic anisotropy on the 
plastic flow properties of the plate. 

Mechanical anisotropy has been recognized as the 
probable cause for the variation with direction of the 
fracturing characteristics of some wrought ,metal 
products.*~? In processing the metal, certain phases, 
inclusions, and cavities may be elongated in the direction 
of the principal strain and thereby constitute a me- 
chanical fibering or mechanical anisotropy. This condi- 
tion may be retained after annealing and heat-treat- 
ment and may also contribute to the directional proper- 
ties of sheet.§ 

The mechanical properties of forgings are probably 
more dependent upon the mechanical fibering than on 
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the crystallographic fibering. It has generally been 
observed>~’ that the yield strength and tensile strength 
of forged carbon steel billets and bars vary little with the 
(per cent) reduction in area by forging and with the 
direction of the test specimen in relation to the fiber. 
However, the elongation and contraction in area in 
tension tests and the impact strength in notched bar 
impact tests increase with increasing reduction in the 
direction of the fiber and may decrease simultaneously 
in the transverse direction. Hard aluminum- and 
magnesium-alloy forgings’ also have lower transverse 
strength than longitudinal strength because of the 
decrease in contraction in area in the test. 

No data are available for plate materials, but it 
would be expected that, as in forging, the mechanical 
fibering would cause a weakening in one or both of the 
transverse directions, and, for rolling, the normal direc- 
tion would be expected to be the least strong. 

This paper will report on the fracturing characteris- 
tics of tensile specimens taken at various orientations 
in commercial 24ST aluminum-alloy plate. 


MATERIAL AND PROCEDURE 


The distribution of fracturing characteristics was 
made on a commercial hot-rolled 1 '/>-in. 24ST aluminum 


plate. 
All tests were made on 1'/»-in. threaded-end tensile 


specimens as shown in Fig. 1. 


. $ 


CONVENIENT 


The orientations of the specimens cut out of the plate 
are shown on the stereographic projection, Fig. 2. 


Points on this projection represent the intersection of | 


straight lines from the center of the stereographic 
sphere with the surface of the sphere. The points may 
then be considered to represent the longitudinal axes of 
the specimens. 

The positions of the specimens can be more easily 
identified by passing planes normal to the rolled surface 


of the plate at increments of 15° to the rolling direction. 
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SVG PROSECTION SHOWING ORIEN - 
TATION OF THE FEST SPECIMENS. POINTS REPRE - 
SENT THE SPECIMEN AXES. 


Thus, in this system plane A is defined by the rolling 
direction and the normal direction, plane B is defined by 
a line in the rolling plane at 15° from the rolling direc- 
tion and by the normal direction, and so forth, to plane 
G, which is defined by the transverse and normal direc- 
tions. Any specimen can then be located by giving the 
plane designation and the angle of elevation of the 
specimen’ axis from the rolling plane. 

A number was stamped and a line was scribed on the 
end of each specimen blank after it was sawed from the 
plate, so that a record of its original position in the plate 
could be kept. All specimens were made in duplicate, 
the test section in all cases being taken in the middle of 
the plate. Some tests were discarded because of faulty 
technique, and some additional tests were made as 
checks. 

All specimens were reheat-treated after machining in 
the following manner: solution treatment—920°F. for 
45 min.; cold-water quench; aged at least 4 days at 
room temperature. This treatment was carried out to 
standardize the heat-treatment of all the specimens and 
to remove any residual effects of machining. 

The specimens were tested in a Riehle hydraulic test- 
ing machine. The specimens were measured after 
fracture to determine the final area, and from this the 
reduction in area was calculated. 

Fracture stress was defined as the breaking load 
divided by the final area.* © 

Almost all of the specimens fractured on a single 
plane (cleavage type) at an angle to the specimen axis, 
The angle between the normal to this fracture plane 
and the specimen axis was defined as the fracture angle. 
This angle was measured on a Wilder microprojector 
at a magnification of 20X. However, a few of the speci- 
mens which fractured on what seemed to be a series of 


* This fracture stress is actually the average stress at fracture, 
since most of the specimens necked. This caused a nonuniform 
distribution of stress over the cross section. 
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parallel planes could not be measured accurately and 
have not been reported. A special series of tests was 
made to check the general symmetry of the plate— 
namely, that the anisotropy in each quadrant of the 
system of coordinates formed by the rolling direction, 
the transverse direction, and the normal direction was 
the same. Only four quadrants needed to be checked 
since the test sections were all taken from the center 
plane of the plate. These tests were made on speci- 
mens located in the four possible C planes at an angle of 
45° to the rolling plane. ‘ 

The data for these tests are givenin Table 1. It can 
be seen from these data that the values of each property 
measured were equal for all four quadrants, within the 
limits of accuracy and scatter. Therefore, only one 
quadrant as represented in Fig. 2 need be considered. 


TABLE 1 


Properties of Specimens from the Four Possible C Planes 
at an Angle of 45° from the Rolling Plane 


Fracture Reduction Fracture 7 
Stress, in Area, Angle, 
Quadrant Lbs. per Sq.In. - Per Cent Deg. 

1 78,600 14 36 
1 73,500 11 41 
a 76,500 13 39 
2 73,100 10 37 
3 73,600 11 
3 71,200 11 35 

4 77,900 12 39 
4 78,100 14 42 


RESULTS AND DISCUSSION 


The fracture characteristics of the aluminum-alloy 
plate were found according to Figs. 3 and 4 to be 
highly directional. While data for sheet and plate for 
directions other than within the rolling plane have not 
been previously determined, it was to be expected from 
investigations on flat-forged products’ that normal 
specimens would be considerably less strong and less 
ductile than those within the rolling plane. 


The data in Figs. 3 and 4 confirm these expectations. 
However, they reveal a rather complicated dependence 
of the fracture characteristics upon the direction of 
testing. These should be correlated in some way with 
the type of break, which was also found to be direc- 
tional. Therefore, the relation between the fracture 
appearance and the orientation may be discussed first. 


The fracture surface of all specimens consisted nearly 
of a single plane.t This plane was generally inclined 
to the specimen axis, the angle between the testing direc- 
tion and the normal to the fracture surface varying be- 
tween approximately 38° and zero, depending upon 
the orientation, as shown in Fig. 5. Furthermore, all 
these fracture plane normals were found to be located 
in the plane defined by the specimen axis and the plate 
normal, as can be seen in Fig. 6. This picture illus- 
trates the fractures for different specimens in the A 


+ As contrasted with cup and cone fractures. 
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plane, which plane contains the rolling direction and 
the normal direction. 
’ The breaks observed were of two distinctly different 
types. In all instances where the fracture angle, in 
Figs. 5 and 6, was distinctly less than 38°, the fracture 
plane was practically parallel to the rolling plane—i.e., 
the plane of weakness for this particular product. 
Consequently, this fracture angle was equal to the angle 
between the specimen axis and the plate normal. How- 
ever, if the specimen axis was inclined more than 38° 
from the plate normal, the fracture angle was found to 
be practically constant—namely, 38 = 5°. 

A macroetched section, Fig. 7, taken perpendicular 
to the rolled surface, shows the presence of a fibrous 
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Fic. 7. Macro-etched section of 1'/s-in. 24ST plate taken- 


perpendicular to the rolled surface. Etched alternately in 5 per 
cent NaOH and Tucker’s etch (2X). 


structure that could cause the presence of a weak plane 
parallel to the rolling plane. 

The values of fracture stress, in Fig. 3, show that there 
is a similar trend in all the planes A to G, which were 
defined as planes containing the plate normal and a 
direction in the rolling plane at an angle to the rolling 
direction. In all the planes a maximum was observed 
at or near the rolling plane. The fracture stress then 
decreased as the orientation of the specimens approaches 
the normal direction. At an angle of 52° from the roll- 
ing plane (38° from the plate normal), the fracture 
stress seemed to level off at a minimum value equal to 
the value for the normal direction. In the rolling 
plane there was a definite trend from a maximum of 
90,000 Ibs. per sq.in. in the rolling direction to a mini- 
mum of 80,000 Ibs. per sq.in. in the transverse direction. 
In general, the fracture stress was constant for orienta- 
tions near the normal direction and increased as the 
orientation moved to the rolling plane, the increase 
being greater in the plane containing the rolling direc- 
tion than in the plane containing the transverse direc- 
tion. 

The ductility (reduction in area), Fig. 5, followed the 
same trend as the fracture stress, except for the rolling 
plane. The value of ductility in this plane remained 
almost constant from the rolling direction to about 
60°, where it changed slope sharply to a minimum at 
90° (transverse direction). This change in behavior 
may be attributed in some way to the crystallographic 
orientation. 

The curves for fracture stress, ductility, and fracture 
angle show breaks in the curves for all planes A to G, at 


. 


1948 


approximately 52°. This would indicate that with the 
presence of a weak plane the process of fracture results 
in competition between planes at a constant angle, 
determined in some manner by the mechanical aspects 
of fracturing, and the plane of weakness. However, if 
the normal stress acting on the fracture plane should be 
considered as the criteriorf of fracture, the fracture 
stress would have to be constant for orientation angles 
from 0° to 52° and would decrease in proportion to a 
secant squared function of the angle from 52° to 90°. 
Actually, the fracture stress and ductility values deviate 
considerably from this behavior. The deviation can- 
not be considered to be due to crystallographic aniso- 
tropy because the effect is three ‘times as large as the 
effect on the yield stress’. 

Therefore, the changes in fracture stress and ductility 
cannot be explained because of the inadequacy of our 
knowledge of the initiation and propagation of fracture. 


CONCLUSIONS 


(1) The fracture angles of the tensile specimens 
indicated the presence of a plane of weakness parallel 
to the rolling plane. A macrosection showed the de- 
velopment of.a fibrous structure that could account for 
the presence of such a plane of weakness. 

(2) The fracture stress and ductility exhibited a 
pronounced directionality. However, this direction- 
ality cannot be correlated with the conception of a 
plane of weakness using any simple fracture law. 
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The Theory of the Impact Tube at Low 
Pressures 


P. L. CHAMBRE? anv S. A. SCHAAF? 
University of California at Berkeley 


ABSTRACT 


A theoretical analysis has been made for an impact tube of the 
relation between free-stream Mach Number and the impact and 
free-stream pressures and densities for extremely low pressures. 
It is shown that the results differ appreciably from the corre- 
sponding continuum relations. 


NOMENCLATURE 
b =1+s 
c = +/2RT, most probable molecular speed 
erf(x) = the error function of x, erf(x) = (2/+/ x) S *e~"dy 
L = length of tube 
M = Mach Number of free stream 
N = molecular number density 
p = pressure 
R = gas constant 
r = radius of tube 


s = U/c,, molecular speed ratio 
T = temperature of the gas 
x 


macroscopic velocity of stream 

coordinate along axis of tube 

W(L/r) = probability of a reservoir molecule escaping through 
a tube of length to radius ratio L/r 

Wrr(X) = probability density that a molecule leaving a ring 
of area 24 r dx will directly fall on another ring 
of area 27 r dx a distance x away without hitting 
the wall 

Wwr(x) = probability that a molecule that leaves a ring of 
of area 27 r dx will go through a cross section of 
the tube at x distance away from the ring with- 
out hitting the wall 

Wer(x) = probability density that a molecule on leaving a 

cross section of the tube will fall directly on a 

ring of area 27 r dx at a distance x away from the 

cross section 
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Wee(x) = probability that a molecule leaving a cross section 
of the tube will directly fall on another cross 
section at a distance x away 


w(x) = the exit probability 

p = density 

Y = ratio of specific heats, taken as 1.4, 
Subscripts 


1 refers to free stream 
2 refers to reservoir 

IS refers to isentropic 
Ra refers to Rayleigh 


(1) INTRODUCTION 


eg IS THE PURPOSE of this paper to provide a theoret- 
ical basis for the interpretation of impact tube 
measurements taken in a high-velocity gas stream of 
such low pressure that the mean free path is larger than 
the impact tube diameter. The familiar Rayleigh for- 
mula, which relates the free-stream Mach Number 
and the ratio of impact to static pressure, is based on 
continuum theory and, hence, becomes invalid under 
these conditions. 

A formula is obtained which relates the free-stream 
Mach Number and the ratio of the product of pressure 
and density in the free stream to that in a reservoir (i.e., 
the pressure-measuring device) located at the end of 
the impact tube. The formula contains a parameter 
dealing with impact tube geometry in contrast to the 
Rayleigh formula, which does not depend on this fac- 
tor. The tube is supposedly placed in a gas stream 
longitudinal axis parallel to the macroscopic velocity 
and one end connected to a reservoir in which the 
pressure-measuring instrument is located. The de- 
rived expression indicates that large errors may be made 
by using the Rayleigh formula for the determination 
of the Mach Number. 

While the chief aim of the theoretical development 
has been to obtain an impact tube formula applicable 
to measurements at low pressures, the method may 
easily be generalized to other free molecular flows of 


interest. 


(2) RESULTS 


It is assumed that (1) all collisions between gas molecules may be neglected (i.e., free molecular flow); (2) gas 
molecules are reflected diffusely from the inside surface of the impact tube with an average velocity characteristic 
of the surface temperature (which is assumed constant); (3) gas’molecules are emitted from the surface at the 
same rate as they are incident (in other words all outgassing phenomena are neglected); (4) a steady state has 
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been attained; (5) the gas in the free stream appears to an observer, traveling with the macroscopic velocity of 

the stream, to be in Maxwellian equilibrium (the gas in the pressure-measuring device is also assumed to be jn 

Maxwellian equilibrium); (6) there is no disturbance in the gas stream (such as a shock for example) due to the 

presence of the impact tube (this appears to be a plausible assumption since the mean free path is extremely large), 
Under these conditions it is shown that 


W(L/r) 
where the discharge fraction W(L/r) is a numerically calculated function, the graph of which is contained in Fig, |, 
The relation indicated in Eq. (1) is presented graphically in Fig. 2 for four different values of the parameter L/y, 


As a means of comparison with continuum theory, the Rayleigh formula! has been expressed in terms of the ratio 
of the product of pressure and density in the reservoir to that in the free stream. The result is 


+ -1) 
P2p2 | 27 y 


and has been included in Fig. 2 as indicated. Eq. (2) assumes a normal shock wave between the free stream and the 
impact tube, followed by isentropic deceleration to zero velocity. Since the shock zone is believed to become 
thicker with increasing mean free path, a modified Rayleigh formula, which omits the shock wave entirely (i.e., 
considers only isentropic deceleration), might also be considered of interest at low pressures. The corresponding 
relation is 


(pope/Pipr)rs = {1 + [(y — 1)/2] + - @) 


and has also been included in Fig. 2. It should be kept in mind that the pressure in the system be everywhere low 
enough so that the assumption of free molecular flow is valid. 

The present theory treats only the equilibrium state. The question of the time lag in the response of such a 
system to a change in pressure has been investigated .for similar free molecular flow conditions.? In particular, 
the advantage to be gained by increasing the impact tube length so as to obtain magnification of the pressure 
density ratio as indicated in Fig. 2 must be balanced against the corresponding increase in the response time lag of 
the instrument. 


(3) THEORETICAL DERIVATION 


In order for a reservoir molecule to escape down the 
impact tube into the free stream, it must either go 
- <~ directly down the tube without hitting the side walls or 

on must first hit at some point a distance x from the reser- 

— voir and then wander out by means of multiple reflec- 

tions from the walls. Hence, the probability of escape 

for a reservoir molecule that has entered the impact 
“ tube is given by:* 


L 
W(L/t) = weL) + dx (A) 


005 where the exit probability w(x) must satisfy the equa- 
tion ; 


L 
on f — de’ (5) 


.00! 

e.! 5 1,0 5.0 1.0 50 100 
The probabilities w.., Wer, Wre, and w,, are obtained 
from simple geometric considerations under the assump- 
DISCHARGE FRACTION wth) AS tion of completely diffuse reflection from the internal 


onto corresponding parts of a circumscribed sphere and 
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calculating area ratios. The results are :* 


A FUNCTION OF TUBE GEOMETRY . surfaces of the tube by projecting the areas in each case — 
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(6) 


Substituting the expressions (6) into Eq. (5), we obtain 
an integral equation of the Fredholm type for the exit 
probability w(x). For a fixed value of L/r, this equa- 
tion may be solved approximately for w(x) in the form 
of a linear function of x. The discharge fraction 
W(L/r) is then calculated from Eq. (4). The result is 
presented graphically in Fig. 1. 

The total number of reservoir molecules which enter 
the impact tube in unit time is given by‘ 


Nowr?V RT2/ (7) 


Hence, the total number of reservoir molecules which 
escape into the free stream in unit time is given by 


Nowr?W(L/r)V RT 2/20 (8) 


In a similar fashion, an expression may be obtained 
for the total number of free-stream molecules which 
enter the reservoir in unit time. This is 


Mer + {1 + x 


or) (L/br)? + 2 


L 
| 
(9) 
For equilibrium the expressions in Eqs. (8) and (9) 
must be equal. Hence, with some algebraic reduction 
and use of the perfect gas law, we obtain Eq. (1). 
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HERMAN F. MICHIELSEN 


The Behavior of Thin Cylindrical Shells After 
Buckling Under Axial Compression 


Fokker Aircraft Factories, Amsterdam, Holland 


SUMMARY 


The fundamental investigations of von Karman and Tsien 
on the buckling of cylindrical shells under axial compression 
are continued. The energy expression is simplified and mini- 
mized with respect to the axial and circumferential wave-length 
parameters. Solution of the equations obtained yields curves 


of the reduced average stress and of the wave dimensions plotted. 


against the reduced average strain. They illustrate the behavior 
of the cylinder during the buckling process. The minimum 
buckling stress is found to be 0.195E(t/R). 


INTRODUCTION 


I 1941, von Karman and Tsien! discussed the buck- 
ling of cylindrical shells. By introducing a diamond 
pattern and nonlinear large deflections, they succeeded 
in proving that, after the theoretical critical load is 
exceeded, decreased loads are required to keep the shell 
in buckled equilibrium. A number of curves corre- 
sponding to different values of parameters which char- 
acterize wave length and width show this decrease 
qualitatively but do not give quantitative results re- 
garding the real load when the dimensions of the 
wave can change freely. It is suggested that the 
lowest energy level determines the values of the param- 
eters, but no exact calculations about that level are 
given. 

One is tempted to consider the envelope of curves of 
oR/Et plotted against eR/t as the actual relation be- 
tween these quantities when no restrictions with re- 
spect to the dimensions of the wave are prescribed. 
Computations show, however, that at small values of 
wand 7 the value of cR/Et becomes negative (tension). 
This result is unlikely. 

At the Sixth International Congress for Applied 
Mechanics Leggett? discussed this matter and made the 
suggestion that the actual curve ¢R/Et vs. eR/t could 
be found by determining the points at which the areas 
subtended by two consecutive curves are equal. In 
Fig. 1, curves corresponding to » = 1.0 and to n = 0.20 
and 0.225, respectively, as computed in reference 1, 
are shown. At point A the areas subtended by the two 
curves are equal and thus A would be a point of the 
real curve if « = 1.0 could be considered as given. 
This, however, is not the case, and when 7 is 
eliminated the same procedure has to be repeated for 
several values of » to eliminate the second parameter 
also. 
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Although the suggestion is correct, * a simpler solution 
seems to be the differentiation of the energy quantities 
with respect to 7 and yw. The differentiation with re- 
spect to f; and f, was performed in reference 1. The 
subject of this paper is, first, a simplification of the 
derivations given in reference 1 and, secondly, differ- 
entiations with respect to 7 and wu. The latter differ- 
entiations were made possible through the simplifica- 
tions in the original derivations. 

Computations showed that the present procedure 
leads to the same result as the ‘‘equal surface” solution 
of reference 2. 

The investigations yield a single relation between 
oR/Et and eR/t. At any value of eR/t the values of 
£1, 2, n, and uw are uniquely determined by the condi- 
tion that the sum of the internal and external energy 
be insensible to small variations in these parameters. 
In Fig. 2 a few other quantities are plotted against 
un’, which seems to be the most suitable basic parameter; 
in Fig. 3 the same quantities are plotted against «R/t, 


* After this paper was submitted the author received a copy of 
“The Behaviour of a Cylindrical Shell Under Axial Compression 
When the Buckling Load Has Been Exceeded,” by D. M. A. 
Leggett and R. P. N. Jones, A.R.C. Reports & Memoranda No. 
2190, dated August, 1942, but just distributed. The results of 
the two investigations are in excellent agreement, but the mathe- 
matical development presented here is only indicated in the 
British paper. 
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Eq. (30) of reference 1 can then be written in the 
following form: 


Wi + W. — Ws; 
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THIN CYLINDRICAL SHELLS AFTER BUCKLING 


Fic. 3. Average compressive stress and wave form, plotted 
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The following parameters, comparable to those de- 
fined in Eqs. (31) of reference 1, can now be introduced: 


R 
o= 


£2 
_ = 4 
n R (4) 


Differentiation of Eq. (3) with respect to g; and ge 
and setting the derivatives equal to zero leads to: 


oR | 16u4 
2 
ut (1+n)? |, 
2 
and 


oR 
= { G+ 20+ gut 


Lad tu 


[4+ 
2(1 + 6(1 — 


With the new symbols 


n?2(1 + (6) 


py?) 


2 


= (7) 
Eqs. (5) and (6) can be written as 
oR 1 (1 + 
16+4 4(——— 4 
Sw? 
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1 1 5 
t 
1+ gu? 
Ky = 4+ (1/2As) — Buy 


From Eggs. (8) and (9) the factor yy/8w? can be elim- 


inated, resulting in and 


Au = 2)2 
Et 1 3 B = 16 
x x, 1+ gn? gtu (12) 
where Dy = 4(1 + p*)/(1 + 


As it seems that in the early stages of buckling u will be extremely small, the accuracy of the numerical com- 
putations will be favorably influenced by eliminating the factor yu‘ from the denominator in Eqs. (8) and (9), 
The terms 1/A, are now to be eliminated: 


Du- 1 


Sw? 4+ (1/2Ay) — [Bu — (Cu/Au) ly 


4+ (1/24) = [Ba — 
Because of Eqs. (10), (11), and (13), one has 


oR 1 Cu-—1 7 
(s- Au ) y+ Bart] + 


Du Cu — 1 
Ax 7 { 4+ (1/2Au) — [Bu — (Cu/Ax)]y 


Use of Eqs. (12) results in 
(Cu — 1)/Au = wX(1 + w?)*/2u4 = (1 + w*)?/2 = Ey (15) 


Hence, the final expression for ¢R/Et is: 


oR 1 


12(1 — v*)y 
(1/2) + 4Au — (AuBs — 
The shortening of the cylinder can be derived from Eq. (23) of reference 1: 
«eR oR 1R 


np 3 2 1 1 


Cuy (16) 


Application of Eqs. (1), (4), and (7) yields 


Fre 


or 


. 
at 
in. 
ty 
Sat 
of 
if 
oR 
Et 
at 
4 
‘ 
+ 
d > l 
4 


(11) 


(12) 


1 (9). 


(13) 


(14) 


(15) 


16) 
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R_ oR (1 _ oR 
(nf)? (| + 2) + (1 +0 (1+ 


From Eqs. (13) and (15) one obtains 


1201 [5 — (8+ Bu — Ex)y + (Bs 


(1/2) + 4A, — — Cay 


Eqs. (16) and (17) represent, in a closed parametric form, the relationship between ¢R/Et and ¢R/t, as shown 
in Figs. 4 and 5 of reference 1. They are an equivalent but simplified version of Eqs. (35), (34), (33), and (39) of 
reference 1, and the cubic corresponding to Eq. (33) has been avoided. Numerical computations indeed yield the 


same results. 


DIFFERENTIATION WITH RESPECT TO 7 AND yp 


As was mentioned in the introduction, the final equilibrium positions of the shell can be found by differentiation 
of Eq. (3) with respect to m and yu, as well as to g; and go, and by setting the derivatives equal to zero. 
Derivation with respect to m leads to: 


1 2 2\2 4),,)2 


From the manipulation Eq. (18) — Eq. (5) — Sw? X Eq. (6), it follows: 


4u4 1 (1 + 


Introduction of the symbols defined in Eqs. (7) and (12) and division by yy yield: 


Solution for w* gives 


Dz 1 1 1 


(21) 
From Eqs. (21) and (13) it follows: 


Eq. (22) can be solved for 1/12(1 — v?)y?: 


1 
1 Cu 


1 


or 


d 
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Ax Au 1 1 
lise Al 12(1 — — + + 1) € -+- ans 


[Du(A uBu Cu) + Auwl(By An) + + 


[5 + 444) + AB. + 44 (5 4 =0 


When y is determined for a given value of « and several values of y, Eqs. (10) and (11) are most suitable for the 
computation of cR/Et. They can be written: 


Et (120 — 
R Diu 1 
bi 24, oA, Boy) | 
ie A i. in which, in accordance with Eqs. (13) and (15), 
[A u(Du — 1)/12(4 — — [5 — (8 + Bu — Eu)y + (Bu — (25) 
a3 (1/2) + (A uBu Cu)y 
The two Eqs. (24) can be used to check the computation. A further check is found in Eq. (21): 
1201 — * — | (26) 
Finally the strain is found from Eq. (17): 
eR/t = (oR/Et) + (1 + w*) + (7/8?) ] (27) 


Eqs. (23) to (27) define the relation between ¢R/Et and ¢R/t at a given value of uy. Now u must be eliminated by 
differentiation of Eq. (3) with respect to yu: 


+ 80%) = | fort oth - 


4 
Et (1+)? (1+ gu)®  (g + (1 + 2’)? 


The difference, Eq. (18) — Eq. (28), is: 


+[ 32u° Je 


0= + 


Substitution of Eqs. (7) yields: 


Livy | | 4u® 


The final form of this condition can be written as: 


2(1 — v®)y? E + ¥ yu? l Sw? 1+ 


16 + 36 + 24+ 2—y* (29 


4 4 
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(25) 
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The solution of Eq. (29) seems to be possible only by 
trial and error. The procedure to be followed is: 

Take a certain value of yu, or rather yu”, and several 
values of y. Eq. (23) then yields, for each value of y, 
the corresponding value of 1/12(1 — »v?)y*, and from 
Eq. (13) the values of y/8w? may be computed. By 
substituting y, 1/12(1 — »v*)y?, and y/8w? in Eq. (29), 
one can determine the value of y that satisfies the 
equation. Substitution of the final values of y, (1/12) X 
(1 — v?)~*, y/8w*, and y in one of Eqs. (24) and in Eq. 
(27) yields oR/Et and ¢«R/t. One point of the curve 
oR/Et vs. eR/t is found in sucha manner. The compu- 
tations have to be repeated for several values of y’. 

Naturally at every point of the curve the corre- 
sponding values of w = V ¥/8(y/8w*) and of » and ¢ 
according to Eq. (7) can easily be determined. Finally, 
Eqs. (4) yield the values gi, g2, and m as functions of 
t/R, determining uniquely the shape of the buckle at 
the stage of loading concerned. 


RESULTS AND CONCLUSIONS 


Computations were carried out in the range 0.01 < yp? 
< 6, and their results are plotted in Figs. 2and 3. Itis 
seen that the load curve has aminimum value (¢R/ Et), in. 
= 0.195 when p? = 0.16. It is remarkable that a maxi- 
mum value of (¢R/Et) mar 0.77 is reached at = 


2.4. It is, however, questionable whether this decreas- 
ing value of ¢R/Et with further increasing «R/t is a 
physical reality. It appears likely that the wave 
pattern assumed at the beginning of this investigation 
is not accurate enough when the loads greatly exceed 
the buckling load, 

In Figs. 2 and 3 the abscissa is u? and «R/t, respec- 
tively, and the ordinates are oR/Et, eR/t, w = go/gu, 
¢ = g:R/t = w,/t (where w, is the half depth of the 
wave), t/Wn = b/WRt (where b is the half width of 
the wave), and Vv nu? = a/V Rt (where a is the half 
length of the wave). In the region most important for 
practical purposes—namely, when 0.1 < yw? < 2— 
there seems to be little variation in w, w,, and’. Only 
the wave length a decreases gradually. The irregular 
shape of the curves is rather unexpected. A remarkable 
but also questionable result is the decrease of ¢ when 
pw? > 0.8. 

Finally, the stresses in the buckled shell can be cal- 
culated from Eq. (20) of reference 1. With the nota- 
tion of Eq. (1) of this paper, Eq. (18) of reference 1 
can be written: 


A= 1/og,°n? — B= 
C = 4gigon” — gi 


D = G = 2gigen H = 16g2*n* 


The stresses in the median surface are then found from Eqs. (20) of reference 1: 


o,R Rp® 2ny gi mx ny 2gig2n” 3mx ny 
(4¥ — 1)o mx Bmx my 
4w? 2 2 R 
g a 
2w m. ony 4w? 2mx 2ny ) 


The bending stresses can be determined from the curvature of the buckles: 

— [Et/2(1 — + v(0°w/dy?) 
— [Et/2(1 — v*)][(O°w/dy*) + v(0%w/Ox*) 
= —(Et/1 + v)(0°w/dxdy) 


ll 


= 

| 


q 
q 
q 
- 


Hence: 


ng 
Et 2(1 — v?) 
unt 


mx 
= 2(1 — v) sin — si 
Ei ( vy) sin R sin 


Egs. (30) and (31) have not yet been evaluated sys- 
tematically. Therefore no curves of the maximum 
stresses can be given, but a preliminary investigation 
led to the following estimate of the stresses: ; 


The maximum compressive stress —o, in the median 


surface is from 2.8 to 3.7 times the average actual com- 


pressive stress, provided 0.2 < < 1.0. It occurs at 
the point, where mx/R = 0° and ny/R = 81° + 82°. 
Tensile stresses and stresses in the circumferential 
direction do not seem to be of importance. 


The maximum bending stress o,,) in the circumfer- 
ential direction is from 1.2 Et/R to 1.6 Et/R at the 
point where mx/R = ny/R = O (maximum inward 
deflection) provided 0.2 < yw? < 1.0. The maximum 
bending stress o,, in the axial direction is about 2.7 
times the average actual compressive stress at the same 


point. 


v mx ny 2mx v 2ny | | 
= ——— ](14+— — cos — cos 
Ei — + R + 4w cos R cos | 


mx ny 2ny 2mx 
+ cos = cos + 4w cos R + 4wy?y cos | 
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R 


(31) 


| 
R 


As the maximum median and the maximum bending 
stress do not coincide, it seems justified to assume that 
the maximum stress is about 5 times the average actual 
compressive stress. This means that, if the latter does 
not exceed 0.2 o.2,* no permanent deformation will 
occur during buckling. 

A systematical investigation of the total stresses in 
the buckled shell, based on Eqs. (30) and (31), might 
yield important further information. 
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Spectrum of Locally Isotropic Turbulence’ 


LESLIE S. G. KOVASZNAY?¢ 
The Johns Hopkins University 


SUMMARY 


The concept of locally isotropic turbulence is used in Kolmo- 
goroff’s sense. Kolmogoroff’s second similarity hypothesis is not 
used, but other hypotheses are suggested for the similarity of 
energy transfer in the energy spectrum of turbulence. The spec- 
trum thus obtained is identical with the —5/3 power law in the 
nonviscous range but falls off rapidly in the viscous high-frequency 
end. The correlation coefficients agree with Kolmogoroff’s pre- 
dictions for both small r and for large r and give explicit results in 
the intermediate range. The only free constant in the present 
theory is Ry, a certain Reynolds Number that corresponds to the 
smallest eddies. Extensive comparison with available energy 
spectrum and correlation measurements have been made. 


INTRODUCTION 


THEORIES describe turbulence in 
terms of various quantities, such as correlation 
functions of various kinds, energy spectra, and prob- 
ability distribution functions. Those quantities are 
interrelated, some by mathematical and some by only 
physical links. 

Taylor! first showed the relationship between a cor- 
relation function and a spectral function. Since then 
there has been some disinclination to base investiga- 
tions on the energy spectrum, since there exists the pos- 
sibility of using a Fourier transform of the correlation 
function. 

From a purely pragmatic point of view there was the 
experimental difficulty of measuring turbulent energy 
spectra; it was easier to study the turbulence through 
quantities like the intensity or correlation function, 
which already had well-established experimental tech- 
niques. More careful investigation shows that the 
measurement of the spectrum during a finite time in- 
terval has a fundamental uncertainty.” 

During and after the war a new approach emerged 
from the work of various authors.*~* That is the idea 
of an ‘“‘eddy cascade,’ with some similarity concepts 
between the eddy stages. The physical picture is 
visualized as a multistage sequence of eddy sizes or fre- 
quencies where the energy is passed from the lower fre- 
quency stages or frequencies by means of the Reynolds 
stresses caused by the smaller scale turbulence. If the 
Reynolds Number of the whole turbulence is high 
enough, the number of ‘‘stages’’ is sufficiently great so 
that different stages may be ‘‘statistically decoupled” 

Received August 2, 1948. 
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from each other. The viscosity plays a role only in the 
smallest eddy stages—i.e., the highest frequencies. 
Assuming dynamic similarity between the stages it is 
possible to derive the energy spectrum (—5/3 power 
law) and the correlation coefficient (2/3 power law) 
for the medium frequency or length range, where the 
eddy cascade is fully developed but neither the low 
frequency end (nonisotropic) nor the high frequency 
end (viscous effects) is too near. 

It is generally believed that a ‘“‘minimum eddy size” 
or minimum eddy Reynolds Number exists, but it does 
not seem to be justified to impose a priori this expecta- 
tion to the theory. In a short note’ the basic ideas of 
present theory were outlined. 

We investigate a turbulent velocity field that is 
locally isotropic in Kolmogoroff’s sense*~* in a four- 
dimensional space-time domain G and which satisfies 
isotropy in the ordinary sense as well. 

The space coordinates are x, y, z, and the time is ¢. 
The turbulent velocity components are u, v, w. The 
average quantities are denoted by bars and by virtue of 
isotropy 


7? = 


gives the definition of r.m.s. value. The definition for 
the average of a quantity a(x, y, 2, t) is 


SS SS Ux, y, 2, t) dx dy dz dt 
— G — 
SSS S dx dy dz dt 
G 
With Kolmogoroff’s first similarity hypothesis, this is 


invariant for parallel translation as well as for rotation 
and reflection. 


THE ENERGY TRANSFER FUNCTION 


The fluctuation kinetic energy per unit mass is 
Er = (1/2) + 0 + = (3/2)u”? 


We might represent the velocity field by a three-dimen- 
sional Fourier-Tieltjes integral, but because of isotropy, 
it is sufficient to study one component only. 

With n, the wave number, we define the energy spec- 
trum E(n) in the following way: Let the contribution 
to %? made by waves having wave numbers between 


and n + dn be E(n), so that 


“= = E(n) dn (1) 
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The total turbulent kinetic energy per unit mass of 
fluid is now 


Ex = (3/2) E(n) dn (2) 


The connection between the energy spectrum E(n) 
and the “longitudinal correlation’? was shown by 
Taylor! to be the Fourier Transform, 


R(r) = E’(n’) cos 24 dn’ (3) 
where n’ is the observed wave number and 
Ri(r) = u(x) u(x + 
R, = u(y) u(y + 


and the connection between the longitudinal and trans- 
versal correlation coefficient for incompressible flow was 
established by von Karman and Howarth:" 


Ri(r) = Ri(r) + (r/2)(dRi/dr) (4) 


Since the Navier-Stokes equations are nonlinear, energy 
is fed from any wave number component to a higher one. 

It is convenient to define an energy transport func- 
tion S(n), which is the amount of energy per unit mass 
that passes in unit time from Fourier components that 
are of wave number less than m to components that are 
of higher wave number than n. 

There is another type of energy transport—viz., the 
viscous dissipation. The total viscous energy dissipa- 
tion per unit mass and unit time, €, has contributions 
from all wave number components: 


ou\? 
= 15» (>) = n*E(n) dn (5) 
oy 0 


If we define 
n 
= 607° f (6) 
0 


The contribution to ¢) between frequencies m and n + 
dn 


de(n)/dn *vn*E(n) (7) 


Similarity Hypothesis for the Energy Transfer Function 


' The energy transfer from one Fourier component to 
another is due to the nonlinear ‘‘inertia terms’ in the 
Navier-Stokes differential equations. The equations 
are linear in both the time derivative and viscous terms. 
The nonlinearity occurs in the form of simple quadratic 
terms. If the Fourier components are phased ran- 
domly, the products of any components of different 
wave number average zero. 

If the turbulence consists of statistically decoupled 
stages, two hypotheses may be suggested instead of 
Kolmogoroff’s second similarity hypothesis: 

First, if the turbulence is locally isotropic in Kolmo- 
goroff’s sense, the energy transfer function S(m) depends 


JOURNAL OF THE AERONAUTICAL SCIENCES—DECEMBER, 


1948 


on E(n) and n only. This hypothesis is based on the 
physical concept that the contribution to S() comes 
from a very narrow frequency band in the immediate 
vicinity of n. 

Introducing the ‘“‘Reynolds Number of a Fourier 
component,” 


R(n) = (1/») V E(n)/n (8) 


the only dimensionally correct form for the energy 
transfer is clearly 


S(n) = (9) 


where K is a nondimensional quantity still depending 
on R(n) 
K = K(R) 


Second, if the turbulence is locally isotropic and the 
frequency spectrum covers a great range, K is a uni- 
versal constant. The hypothesis means that the mech- 
anism of energy transfer is similar in various wave 
number ranges and it does not depend on concurrent 
viscosity effects. This independence of K with respect 
to viscosity effects seems reasonable because of the 
linearity of the viscous terms in the equations of motion. 


THE ENERGY EQUATION 


_The concept of locally isotropic turbulence in Kol- 
mogoroff’s sense includes steadiness in the time within 
its space-time domain, G. If there is a continuous 
supply of energy from the low-frequency end of the 
spectrum, the steadiness is rigorously true. In a de- 
caying isotropic turbulence the steadiness applies only 
for a short time, and the low frequency ends suffer 
energy drain that does not fit into the concept of local 
isotropy. 

For the steady case, the difference of the energy trans- 
fer at frequency and n + dn is equal to the viscous 
energy dissipation in the same interval 


S(n) — S(n + dn) = (de/dn) dn (10) 
Substituting from Eqs. (7) and (9), Eq. (10) becomes 


or, if we write 
= Ro (12) 
then, 
dE 5E 8» 
8 


If we neglect viscosity, we get the ‘‘—5/3 power law’: 
if v = 0, dE/dn = —(5/3)(E/n), and E = Ey(n/m)~", 
where Eon;”* is an arbitrary constant. With the sub- 
stitution of z? = E/n, Eq. (13) becomes 


dz 42 
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The solution of Eq. (14) is 
z [C — (v/Ro)n*] (15) 


where C is an arbitrary constant. Introducing a new 
constant % defined by C = vno’*/Ro, we can write Eq. 


(15) as 
or 
E(n) = — (m/mo)*}* (17)* 
where 
Eo = v?no/Ro* (18) 
The relationship between Eo, mo, Ro is 
Ro = (1/) VEy/no (19) 


and this gives the interpretation of the universal con- 
stant Ry as a characteristic Reynolds Number. [See 
Eq. (8).] 

The expression for the Reynolds Number of a Fourier 
component 


R(n) = Ro(n/mo) ~*”" — Ro (20) 


gives the possibility of determining Ry from experi- 
mental data. 


The Dissipation of Energy 


Kolmogoroff describes locally isotropic turbulence in 
terms of only two independent quantities, viz., € the 
total rate of dissipation, and v the kinematic viscosity. 

By dimensional reasoning, he defines a characteristic 
length » and a characteristic velocity u;: 


n= ye, 


(21) 
(22) 


Consequently, Kolmogoroff’s characteristic Reynolds 
Number is unity: 


uxn/v = 1 (23) 

Choosing our characteristic length as the reciprocal 
of the characteristic wave number and our character- 
istic velocity in such a way that the resulting Reyn- 
olds Number would be just Ro, we define 


no = 1/m (24) 


and 


* Present study is effectively dimensional reasoning with 
E(n) and n as basic dimensional quantities. The result of this 
reasoning is Eq. (17) with two free constants EZ» and m. It is 
further assumed that E’(n’) has the same form as E(n) but with 
different constants; therefore it is used indiscriminately in com- 
parison with the experiments. 

Although Eq. (5) is valid both for E(n) and E’(n’), Eq. (7) 
gives the accurate contribution to viscous dissipation only with 
E(n). This fact was pointed out to the author by G. K. 


Batchelor. 


= Rovng (25) 


The two kinds of definition (of characteristic length 
and of characteristic velocity) are connected by the 
dissipation relation: Substituting E(m) from Eq. (17) 
into Eq. (5), we get 


& = 152? (26) 


or 


Bk 1 (27) 
% 

v3} 

Matching Eqs. (21), (22) and (26), (27), we obtain the 
relationship between the present characteristic quan- 
tities and those of Kolmogoroff : 


no = (28) 
uo = (29) 
The energy transfer function S(m) then becomes 
S(n) = e[1 — (n/mo)”}? (30) 
The rate of viscous dissipation is 


de/dn gives the viscous dissipation per wave number. 
It may be of more interest to compute the dissipation 
per ‘‘stage’’ of eddy cascade. 


Let dn = ndé; then, 
de de n | | 
— = 1-—{- 32 
dé in (“) No (32) 
The maximum of the viscous dissipation per wave 
number is at 


n/m = 1/27 = 0.1925 
The maximum viscous dissipation per stage (per oc- 
tave) is at 

n/n) = 1/W/27 = 0.439 
Fig. 1 shows S(n)/eo plotted against the logarithmic 


scale of n/m. Evidently, the viscous effect is confined 
mostly to the frequency interval (m/10) < n < mo, 
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Fig. 2 shows (mo/¢€0)(deo/dn), the dissipation rdte per 
wave number. Fig. 3 shows (n/é)(deo/dn), the dis- 
sipation rate per stage. 


Correlation Cofficient 


The turbulent velocity field can be described sta- 
tistically by either spectra or correlation functions. 
The double correlation in a turbulent velocity field can 
be specified, in general, as a second rank tensor. If 
isotropy in the ordinary sense is assumed, the tensor is 
given by only two scalar functions R,(r) and R2(r). 

Kolmogoroff introduces a different type of double 


correlation. His two characteristic scalar functions are 
Baa = [user + 1, Xp Xx) — xy xe) (33) 
Brn = Xj, Xx +17) — xj, (34) 


where k # 7.. The locally isotropic turbulence that 
Kolmogoroff introduces deals, in effect, with an in- 
finite sequence of eddy cascades, of which he confines 
himself to a small domain. This naturally de-em- 
phasizes the importance of the quantity #*; for turbu- 
lence that is only locally isotropic, it even may not be 
bounded, since the flow may have greater and greater 
eddies, ad infinitum, which have no effect on relations 
derived in the bounded domain of local isotropy. 


The spectrum found in Eq. (17) has t' e character 
that increases indefinitely with decreasing n. Ifn’— © 
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Similarly, the correlation coefficients R; and R, be. 
come indeterminate when r — 0. However, the Kol- 
mogoroff type of correlation functions, By and B,,, stil] 
remain bounded, so for convenience we define 


Bil(rny) = (35) 
B2(rmo) = (36) 


which differs only by numerical constants from Kol- 
mogoroff’s Bya(r/n) and Byn(r/n). 


In incompressible flow from Eq. (4) 


B2(rno) = Bi(rno) + (37) 
Introducing 

n/no = k (38) 

my = p (39) 


we express 8;(p) by Eqs. (3) and (17): 
Bi(p) = — — cos 24kp) dk (40) 
0 


Unfortunately, the integration cannot be performed 
in closed form. Introducing 


h = kp (41) 
and defining 


Ji(p) (1 — cos dh (42) 
0 
Jx(p) = — cos dh (48) 
0 


0 


B, and becomes 
= —Js+Js (45) 
= (4/3)Ji — (2/3)Je (46) 
Since 
sin 1 — cos 2mp (47 
2 (2xp)? 


but J; and J2 are evaluated numerically. 
For large values of p 
Jy = Cp” — (3/2) (48) 
Jz 3 — (49) 


where 


C -f h~** (1 — cos 2h) ch (50) 
0 


M 
C, = 2 (1 — cos 2h) dh (51) 
7 0 


where 
extreiT 


For 
power 


and 


dn 
20 
poly. 
but 
0.01 
(0.02 
0.03 
0.04 
0.05 
a 0.06 
0.08 
ep 0.10 
0.15 
0.20 
0.25 
0.30 
0.35 
0.40 
0.45 
0.50 
0.55 
0.60 
0.65 
0.70 
0.75 
0.80 
0.85 
0.90 
0.95 
1.00 
1.10 
1.15 
1.20 
1.25 
1.30 
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TABLE 1 
— 
2m om bom 
2 +4.935 +9.870 
4 —2.473 —7.422 
6 +1.068 +4.273 
8 —().337 —1.685 
10 +0.079 +0.472 
12 —0.014 —0.098 


where / is chosen so that 1 << M << p. Thus, for 
extremely large values of p, 


C sin?rp 
Bile) = Cp +3,-4- (52) 
p 
2C; 


4 2 
Bo(p) Cp + | 4 (53) 
3 3p 


For small values of p, we can expand 8; and £2 in 
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1 
(1 = k’’)*dk 
0 
—— (55b) 
(3m — 1) (8m + 1)(m + 1) 


(—1)™+? (2x)™ 


(2m)! (3m — 1)(3m + 1)(m + 1) 


(56) 


The series is alternating and convergent for any value 
of p, since 


Lt 


From Eq. (37) it is evident that 
bom = (1 + mM) dom (57) 


The functions 8; and #2 satisfy the predictions of 
Kolmogoroff for both p<< land p >> 1: 


power series. Both of them are even functions: Ifp >>, 
Bi = Cip 
Bi(p) = (54) 4 
m=1 2/5 Bo 
and 1 
If p<<l, 
oat 2m 55 
B2(p) Pimp (55) By = as be = 2ae 
2m 1 Be 
dp?” 0 Bi 
(55a) The constants for the asymptotic expansion are: 
but C, = 6.844, C, = 0.398 
TABLE 2 
p By Be p Bi Be p Bi Be p Bi Be 
0.01 0.0005 G.0610 | 1.35 4.704 7.356 3.9 13.12 18.72 | 13.0 33.91 46.50 
0.02 0.0020 0.0039 | 1.40 4.900 7.619 4.0 13.40 19.10 | 14.0 35.82 49.05 
0.03 0.0044 0.0089 1.45 5.093 7.881 4.1 13.69 19.48 i; 15.0 37 .69 51.54 
0.04 0.0079 0.0158 | 1.50 5.264 8.142 4.2 13.97 19.86 16.0 39.52 53.98 
0.05 0.0123 0.0246 | 1.55 5.474 8.402 4.3 14.25 20.23 | 17.0 41.31 56.37 
0.06 0.0177 0.0354 | 1.60 5.662 8.661 4.4 14.53 20.60 18.0 43.06 58.71 
0.08 0.0364 0.0629 | 1.65 5.850 8.920 4.5 14.80 20.97 19.0 44.79 61.01 
— 1.70 6.035 9.176 | 4.6 15.08 21.33 - 
0.10 0.0491 0.0980 | 1.75 6.220 9.431 | 4.7 15.35 21.70 20.0 46.48 63.27 
9.15 0.1098 0.2184 | 1.80 6.402 9.683 | 4.8 15.62 26 | 

, > | 1.95 6.943 10.425 5.0 16.15 22.78 26.0 56.11 76.11 
0.30 0.4248 0.8312 28 0) 59.15 80°17 
0.35 0.5693 1.105 2.00 7.121 10.660 5.2 16.67 23.48 | 399 62.12 84.13 
0.40 0.7304 1.406 5.4 17.19 24.17 39 () 6502 8800 
0.45 0.9062 1.727 | 2.10 7.450 11.185 | 5.6 17.71 24.86 | 67.86 «91.79 
0.50 1.095 2.064 | 2.20 7.816. 11.600 | 5.8 18.22 25.54 36 70°65 95.61 
0.55 1.295 2.412 2.30 8.157 12.050 | 6.0 18.72 26.21 38 0 73°39 99.16 
0.60 1.508 2.765 2.40 8.493 12.500 | 6.2 19.21 26.87 | 
0.65 ee 3.111 2.50 8.827 12.945 | 6.4 19.71 27.53 40.0 76.08 
0.70 1.936 3.474 2.60 9.156 13.390 | 6.6 20.19 28.18 |- ———___— 

0.75 2.157 Se a 9.480 13.830 | 6.8 20.67 28.83 45.0 82.61 
0.80 2.380 4.161 | 2.80 9.801 14.260 | 50.0 88.91 
0.85 2.602 4.493 | 2.90 10.120 14.690 | 7-0 21.15 29.46 55.0 95.00 
0.90 2.824 4.814 | 3.0 10.43 15.11 75 22.33 “31.03. | 80.0 100.91 
0.95 3.044 5.125 | 3.1 10.74 15.53 8:0 23.47 32.57 | 85.0 106 66 
1.00 3.261 5.427 | 3.2 11.05 15.94 | Q's 24 60 34.07 70.0 112.26 
1.05 3.476 5.720 | 3.3 11.35 16.35 75.0 117.73 
1.10 3.688 6.005 3.4 11.65 16.75 9.0 25.70 35.54 80.0 123.08 
1.15 3.897 6.284 | 3.5 11.95 17.15 | 128.32 
1.20 4.103 6.557. | 3.6 12.25 17.54 10.0 27.85 38.41 90.0 133.46 
1.25 4.306 7.226~..1.3.2 12.54 17.94 11.0 29.93 41.19 95.0 138.50 
1.30 4.506 7.002 | 33 12.83 18.33 12.0 31.95 43.88 100.0 143.46 


Kol- 
» Still 
= | 
(38) { 
(39) 
40) | 
med ! 
(41) 
43) 
44) 4 
45) 
46) 
47) 
18) 
50) 
51) 


750 


For the power series, the numerical values are shown in 
Table 1. 

The functions 8; and (2 are tabulated in Table 2. 

The idealized physical picture is an infinite eddy 
cascade with a constant rate of energy transfer between 
the stages in the range of frequencies where viscosity 
effects are negligible. 

in the actual turbulence, there is probably a limit 
with respect to maximum eddy size (or minimum fre- 
quency) related to the dimensions of physical boundar- 
ies, so that the laws found above cannot apply to the 
low end of the spectrum. If a turbulent flow includes a 
substantial intermediate frequency region for which the 
dissipation is negligible and local isotropy still holds 
(this is possible if the overall turbulence Reynolds 
Number is high enough), then no contribution is made 


to the viscous dissipation by those still larger eddies . 


(lower frequencies) that do not obey local isotropy. 
In this case we can express the relationship between \, 
the microscale of turbulence, and m: Introducing 


2 2 
dr* Loo dr* 


and a Reynolds Number R) connected with it, 
Ry = X\u'/v (59) 


From Eqs. (3) and (58), we can deduce Taylor’s now 
familiar expression for the rate of dissipation, ¢ = 


° 
w 


i.e., to lower and lower frequencies. 
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l5vu*A~*. After this expression has been combined 
with Eq. (26), Eq. (59) may be written as 
Ry = (Amo)? (60) 
Introducing 
R 
(61) 
We get 
mo = x/d (62) 
and 
p = x(r/d) (63) 


From Eqs. (35) and (36), we get expressions related to 
those given by Kolmogoroff * 


Ri(r) =1- (uo?/ By (p) (64) 
R2(r) = 1 — (uo?/it?) Bo(p) (65) 


Since = R,/R; [from Eqs. (25) and (60)], 
Eqs. (64) and (65) become 


(66) 


r 1 r 
This shows that the correlation curves plotted against 
r/N are a single family with the single parameter x or R,. 
Fig. 4 shows the single family curves for R; and Fig. 5 
for R2, both plotted against r/X. 

If \ itself is not known directly, we may represent the 
correlation functions in a somewhat different way. We 
can make the distance r dimensionless by combining it 
into a Reynolds Number with the r.m.s. velocity fluctua- ° 
tion, instead of simply dividing by X: 


R, = ru'/v (68) 
and 
R ru 
= —— 69 
(68) 
Then, from Eqs. (63) and (69) we can show that 
p = (70) 
In these variables, R; and R; become 
R(t) = 1 — (71) 
mx 
1 
R.(é) = 1 — hal€/x) (72) 


Figs. 6 and 7 show R,(¢) and R,(é) plotted against &. 
The whole computation of correlation coefficient is 
based on the assumption that the eddy cascade is ex- 
tended indefinitely to greater and greater distances— 
This cannot be 
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true beyond a certain value of r which was estimated 
by Batchelor.'? Thus, the curves obtained here ob- 
viously are comparable to the experiments only at high 
values of the correlation coefficient. 


Comparison of the Theory with Experimental Data of 
Energy Spectrum 

The few published spectrum measurements have 
been made at comparatively low Reynolds Numbers. 
Therefore, only the upper frequency range of the spec- 
trum can be expected to satisfy the requirements of local 
isotropy. 

The simplest method of comparison is in terms of the 
Reynolds Number of the Fourier components, accord- 
ing to Eq. (20). If we write n~‘* = pand mm” 
Po, then 


R(p) = (Ro/po)p — Ro (73) 


This means R plotted against p must be a straight line 
intersecting the p axis at pp and the R axis at —R». 
The measured spectra of Simmons™ are used for com- 
parison with the present theory. We denote 


f 


frequency, cycles per sec. 


U, = mean speed 

u’ =r.m.s. turbulent fluctuations 

F(f) = Simmons’ energy spectrum function, which 
is normalized by 


= 1 
0 


(74) 


We transform 
n = f/Uo 
and 
E(n) = Uu"F(f) 


Figs. 8, 9, 10, and 11 show Simmons’ measurements R 
plotted against p = n~‘. The agreement is good and 
the points lie on a straight line. It is assumed that 
v = 0.13 cm. sec.~*. The value of Ro is found to be 
11.5. 

Naturally the measurements are not too accurate 
in the upper end of the frequency range where the noise 
effect is great. Using the valties of po, we get mo and the 
corresponding maximum frequency fo = Uomo. 

The turbulence spectrum was measured 120 in. down- 
stream of a 3-in. mesh length grid, and the data are 
given in Table 3. 


(75) 


TABLE 3 
Us no to 
Ft. per sec. Cm.~! Cycles per sec. 
35 3.70 3,975 
30 3.40 3,110 
25 2.95 2,250 
20 2.50 1,525 
120 T 
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Using the values obtained for mp) and plotting E(n)/ 
Ko against n/no, we find good agreement between the 
experimental points and the theoretical curve, Fig. 12. 
The deviation at lower frequencies is attributed to the 
fact that the local isotropy probably is not valid for a 
domain large enough to include the lowest fre- 
quencies; this part of the spectrum is governed by 
another physical picture, .t covered by the present 
theory. 


Comparison of the Theory with Correlation Measurements 


The family of correlation curves (Figs. 4, 5, 6, and 
7) may be checked with experiments if the conditions 
of turbulence correspond to the postulated local iso- 
tropy. 

Although in the Fourier transform process from 
spectrum to correlation curve all regions of the spec- 
trum affect all regions of the correlation function, this 
influence is not of equal ‘‘weight.”’ 

The total energy dissipation ¢) and the value of \ 
both get their principal contributions from the upper 
end of the spectrum. The power series expansion of 
the correlation coefficient around the origin (ry = 0) 
shows [Eq. (55a)] that the contribution to the higher 
power terms is relatively small from the low-frequency 
end. 
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The turbulence behind regular grids seems to have ag 
rather characteristic spectrum shape in the low-fre. 
quency end. Dryden,'* using experimental  eyj- 
dence, assumed the correlation function R; to be 
exponential. This assumption, transformed into 
spectral representation, results in a simple form of the 
spectrum: 


E'(n’) = 40°L,/[1 + (2xn’L,)?] (76) 
where 
L,= Ri(r) dr (77) 
0 
or 
L, = E(0)/4i (78) 


Naturally, if we were to construct a “‘theoretical’’ cor- 


relation curve by fitting the present theory for the 


high-frequency (small r) range with the empirical curve 
[Eq. (76) ] in the low-frequency range, this would check 
experiment better than the high-frequency theory 
alone. This kind of treatment has been avoided since 
it does not contribute to the validity of the present 
theory. 

Corrsin has made available some unpublished cor- 
relation measurements made at the California Institute 
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of Technology in 1942.'* Fig. 13 represents a set of this 
correlation data, measured behind a 1-in. mesh length 
grid consisting of two perpendicular layers of d = '/,-in. 
rods. The value of Ry = 11.5 was assumed. 
Table 4 gives details. x») is the distance down- 
stream from the grid. No wire length correction has 


been made. 


TABLE + 
xo/M Rx \/M 
18!/. 1.4141 72.3 0.105 
423/, 1.366 67.7 0.167 
67!/» 1.262 57.6 0.206 
91!/> 1.219 53.7 0.231 
115'/, 1.213 53.3 0.267 


We see that the points agree well with the theoretical 
curve computed on the assumption of Ry = 11.5. Fig. 
13 shows the upper part of correlation curves. (R: > 
0.7.) As one example, the xo/1/ = 42*/,, the full cor- 
relation curve is shown in Fig. 14. It is clear that for 
large values of 7, the theory cannot be expected to agree 
with these experimental results. 


CONCLUSION 


A theory is given to describe the spectrum of locally 
isotropic turbulence both in the “nonviscous’’ region 
where Kolmogoroff’s second similarity hypothesis may 
apply and in the viscous region where the dissipation 
has an important role. 

For large eddies compared with dissipative eddies or 
(i.e., low f), the agreements and discrepancies with 
experiment are the same as for Kolmogoroff’s theory. 
For the viscous region the present theory gives explicit 
expressions with only one empirical numerical constant, 
Ry, for the spectrum and correlation functions. The 
agreement with experiments is satisfactory but could 
be greatly improved by some appropriate theory for the 
mechanism of energy distribution at the low-frequency 
end of the spectrum.* 

* Th. von Karman has recently presented (paper read at the 
meeting of the Wes‘ Coast Heat Transfer and Fluid Mechanics 
Institute, June 23, 1948) a generalization of the modern turbu- 
lence theory, valid in the complete ‘‘nonviscous”’ range of the 


energy spectrum. 
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